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1.0  SCOPE 

1.1  Identification 


The  software  described  in  this  document  is  identified  as  the  Polar  Ice  Prediction  System  (PIPS) 
Version  3.0.  PIPS  3.0  is  a  dynamic  sea-ice  model  that  forecasts  conditions  in  all  sea-ice  covered 
areas  in  the  northern  hemisphere  (down  to  30°  north  in  latitude).  It  has  a  horizontal  resolution  of 
approximately  9  km.  The  vertical  resolution  in  the  model  has  been  set  at  45  levels  so  that  Arctic 
shelves,  continental  slopes  and  submarine  ridges  are  accurately  represented.  This  allows  for  17 
levels  in  the  upper  300  m  of  the  water  column  and  a  maximum  layer  thickness  in  the  deep  ocean 
of  300  m.  The  array  size  is  1280x720.  Currently  the  domain  includes  the  Irminger,  Labrador, 
North  and  Baltic  Seas  on  the  Atlantic  side  and  the  Bering  Sea,  Sea  of  Japan  and  the  Sea  of 
Okhotsk  on  the  Pacific  side. 

PIPS  3.0  model  bathymetry  south  of  64°  N  is  derived  from  the  ETOP05  database,  Navy 
Research  Lab  charts  and  Canadian  Hydrographic  Service  charts.  Bathymetry  north  of  64  N 
comes  from  the  2.5  km  resolution  digital  International  Bathymetric  Chart  of  the  Arctic  Ocean 
(IBCAO). 

The  PIPS  3.0  system  uses  the  Los  Alamos  ice  model,  CICE  (version  3.1),  containing  improved 
procedures  for  model  thermodynamics,  physics  parameterizations  and  energy  based  ridging.  It 
has  the  ability  to  predict  multi-category  ice  thickness.  The  CICE  model  is  designed  to  run  on 
massively  parallel  computers  [37,10,11].  The  CICE  model  is  presently  being  coupled  (via  file 
transfer)  to  the  operational,  global  Navy  Coastal  Ocean  Model  (NCOM),  to  predict  ice  thickness, 
concentration,  and  drift  in  the  Arctic  Ocean.  NCOM  is  a  baroclinic,  hydrostatic,  Boussinesq, 
free-surface  ocean  model  that  allows  its  vertical  coordinate  to  consist  of  sigma  coordinates  for 
the  upper  layers  and  z-levels  below  a  user-specified  depth  [5], [34].  NCOM  runs  operationally  at 
NAVOCEANO  at  a  resolution  of  1/8°  globally.  PIPS  3.0  also  forecasts  surface  ocean  current 
and  temperature  in  the  surrounding  seas.  It  is  currently  being  tested  for  success  in  coupling  with 
the  next  generation  global  HYbrid  Coordinate  Ocean  Model  (HYCOM).  For  more  information 
on  HYCOM,  visit  the  website  at  http://hvcom.rsmas.miami.edu/hycom/. 

PIPS  3.0  is  driven  by  heat  fluxes  and  surface  winds  from  the  Navy  Operational  Global 
Atmospheric  Prediction  System  (NOGAPS).  Daily  updates  are  accomplished  through  an 
objective  analysis  of  ice  concentration  data  from  the  Special  Sensor  Microwave/Imager  (SSM/I) 
located  on  the  Defense  Meteorological  Satellite  Program  (DMSP)  satellite. 

There  are  four  primary  components  that  work  together  to  comprise  the  PIPS  3.0  model: 

•  a  thermodynamic  model  that  calculates  snowfall  as  well  as  local  growth  rates  of  snow 
and  ice  due  to  vertical  conductive,  radiative  and  turbulent  fluxes; 

•  an  ice  dynamics  model,  which  predicts  the  velocity  field  of  the  ice  pack  based  on  a 
model  of  the  material  strength  of  the  ice; 

•  a  transport  model  that  depicts  advection  of  the  aerial  concentration,  ice  volumes  and 
other  state  variables; 

•  a  ridging  parameterization  that  transports  ice  among  thickness  categories  based  on 
M  anuscript  approved  October  3,  2008. 
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energetic  balances  and  rates  of  strain. 

1.2  Document  Overview 

The  purpose  of  this  Software  Design  Description  (SDD)  is  to  describe  the  software  design  and 
code  of  the  Polar  Ice  Prediction  System  (PIPS)  Version  3.0.  Because  the  PIPS  3.0  model  is 
largely  based  on  the  CICE  model,  this  document  reflects  the  information  found  in  the  Los 
Alamos  Sea  Ice  (CICE)  Software  User’s  Manual  [1].  The  SDD  includes  the  mathematical 
formulation  and  solution  procedures  for  PIPS  3.0  as  well  as  flow  charts  and  descriptions  of  the 
programs,  modules,  and  subroutines.  This  document,  along  with  a  User’s  Manual  [2]  and  a 
Validation  Test  Report  [3],  forms  a  comprehensive  documentation  package  for  the  PIPS  3.0 
model  system. 

Generally  speaking,  subroutine  names  are  given  in  italic  and  file  names  are  boldface  in  this 
document.  Symbols  used  in  the  code  are  typewritten,  while  corresponding  symbols  in  this 
document  are  in  the  math  font  which  is  similar  to  italic. 
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3.0  PIPS  3.0  SOFTWARE  SUMMARY 


PIPS  3.0  is  written  in  fixed-format  FORTRAN90  and  runs  on  UNIX  host  platforms,  including 
SGI  Origin  3000  (<OS>=  IRIX64  below),  SGI  Altix  (Linux),  IBM  Power4  (AIX)  and  Cray  XI 
(UNICOS).  The  code  is  parallelized  through  grid  decomposition  with  MPI  for  message  passing 
between  processors,  with  four  processors  allocated  to  each  hemisphere.  The  code  has  been 
optimized  for  vector  architectures  and  tested  on  Fujitsu  VPP  5000,  Cray  XI,  and  NEC  platforms. 
At  NAVOCEANO,  PIPS  3.0  is  run  using  32  processors  on  an  IBM  platfonn  with  a 
ConsumableMemory  of  500  mb.  With  respect  to  hardware  resources,  a  one-day  run  of  PIPS  3.0 
at  NAVOCEANO  requires  0.85  Processor  Hrs. 


4.0  PIPS  3.0  SOFTWARE  INVENTORY 

4.1  PIPS  3.0  Components 

A  complete  description  of  each  PIPS  3.0  module  and  subsequent  subroutines,  including  a 
definition,  purpose,  relationship  to  other  modules,  etc.  is  found  in  Section  5.3.  The  PIPS  3.0 
software  is  run  through  a  series  of  modules,  makefiles,  input  files,  and  scripts. 


4.1.1  PIPS  3.0  Modules 

CICE.f,  icealbedo.f,  iceatmo.f,  icecalendar.f,  iceconstants.f,  icecoupling.f, 
ice_diagnostics.f,  ice_domain.f,  ice_dyn_evp.f,  ice_exit.f,  ice_fileunits.f,  ice_flux.f, 
ice  flux  in.f,  ice  grid.f,  ice  history.f,  ice  init.f,  ice  itd.f,  iceitdlinear.f,  icekindsmod.f, 
icemechred.f,  icemodel  size.!',  icempiinternal.f,  ice  ocean.f,  ice  read  write.f,  ice  scaling.f, 
ice_state.f,  ice_thenn_itd.f,  ice_thenn_vertical.f,  ice_timers.f,  ice_transport_mpdata.f, 
icetransportremap.f,  ice  work.f 


4.1.2  PIPS  3.0  Subroutines 

abort_ice.f,  absorbed_solar.f,  add_new_ice.f,  add_new_snow.f,  aggregate. f,  aggregate_area.f, 
albedos. f,  asum  ridging.f,  atmoboundarylayer.f,  bound.f,  bound  aggregate.f,  boundijn.f, 
bound  narr.f,  bound  narr  ne.f,  boundstate.f,  bound  sw.f,  calendar.f,  checkmonotonicity.f, 
checkstate.f,  columnconservationcheck.f,  columnsum.f,  columgrid.f, 

completegetfluxocn.f,  conductivity. f,  conservationcheckvthenno.f,  conservedsums.f, 
construct  fields.f,  departure_points.f,  dumpfile.f,  end  run.f,  evp.f,  evp  finish.f,  evp_prep.f, 
exitcoupler.f,  file_year.f,  fitline.f,  fluxintegrals.f,  freeboard. f,  fromcoupler.f, 
frzmltbottomlateral.f,  getsum.f,  getflux.f,  globalconservation.f,  globalgather.f, 
global  scatter.f,  ice  bcast  char.f,  ice  bcast  iscalar.f,  ice  bcast  logical.f,  ice  bcast  rscalar.f, 
icecouplingsetup.f,  iceglobalrealmimnax.f,  iceopen.f,  iceread.f,  icestrength.f, 
ice_timer_clear(n).f,  ice_timer_print(n).f,  ice_timer_start(n).f,  ice_timer_stop(n).f,  ice_write.f, 
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icewritehist.f,  icecdf.f,  initcalendar.f,  initconstants.f,  initcpl.f,  initdiagnostics.f, 
init  diags.f,  init  evp.f,  initflux.f,  init  flux  atm.f,  init  flux  ocn.f,  initgetflux.f,  initgrid.f, 
inithist.f,  inititd.f,  initmassdiags.f,  initmechred.f,  initremap.f,  initstate.f, 
initthermovars.f,  initthermovertical.f,  init_vertical_profile.f,  inputdata.f,  integer  function 
lenstr(label).f,  interpcoeff.f,  interpcoeffmonthly.f,  interpolatedata.f,  lateralmelt.f, 
limitedgradient.f,  linearitd.f,  loadtracers.f,  localmaxmin.f,  locatetriangles.f, 
make_masks.f,  makemask.f,  merge_fluxes.f,  mixed_layer.f,  mpdata(narrays,phi).f, 
NCAR_bulk_dat.f,  NCAR_files.f,  pipsgrid.f,  popgrid.f,  prepare_forcing.f,  principal_stress.f, 
print  state.f,  read  clim  data.f,  read  data.f,  real  function  ice  global  real  maxval.f,  real  function 
ice  global  real  minval.f,  real  function  ice  global  real  sum.f,  rebin.f,  rectgrid.f,  reduce  area.f, 
restartfile.f,  ridge_ice.f,  ridge_prep.f,  ridge_shift.f,  runtime_diags.f,  scale_fluxes.f, 
scale_hist_fluxes.f,  setup_mpi.f,  shift_ice.f,  sss_clim.f,  sss_sst_restore.f,  sst_ic.f,  stepu.f,  stress. f, 
surfacefluxes.f,  t2ugrid.f,  temperaturechanges.f,  thermoitd.f,  thennovertical.f, 
thickness  changes.f,  timers. f,  Tlatlon.f,  to  coupler.f,  to  tgrid.f,  to  ugrid.f,  transport  mpdata.f, 
transportremap.f,  trianglecoordinates.f,  tridiagsolver.f,  u2tgrid.f,  unloadtracers.f, 
update  fields.f,  update  state  vthermo.f,  zap  small  areas.f 


4. 1. 3  PIPS  3. 0  Input  File  Namelist  Parameters 

The  ice_in.txt  file  is  an  input  file  of  namelist  parameters  used  in  running  the  PIPS  3.0  model  at 
NAVOCEANO.  The  variables  are  described  in  the  order  they  appear  in  the  ice_in  and  with  the 
default  values  and  directory  paths  used  in  execution  at  NAVOCEANO. 

4. 1. 3. 1  Ice  Namelist  (ice  nml) 


Name 

Type/Options 

Description 

Default  Values  /  Directory 
Location 

year  init 

yyyy 

The  initial  year,  if  not  using 
restart 

=  0001 

istepO 

integer 

Initial  time  step  number 

=  0 

dt 

seconds 

Thenno/transport  time  step 
length 

=  2800. 

ndte 

integer 

Number  of  EVP  subcycles 

=  120 

npt 

integer 

Total  number  of  time  steps 
to  take 

=  70 

diagfreq 

integer 

Frequency  of  diagnostic 
output  in  dt 

eg.,  10  is  once  every  10 
time  steps 

=  30 

histfreq 

y,  m,  w,  d,  1 

Write  history  output  once  a 
year,  month,  week,  day,  or 
every  time  step 

=  ’h’ 
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Name 

Type/Options 

Description 

Default  Values  /  Directory 
Location 

dumpfreq 

y,  m,  d 

Write  restart  every 
dumpfreq_n  years,  months, 
days 

=  ’d’ 

dumpfreq_n 

integer 

Frequency  restart  data  is 
written 

=  1 

hist_avg 

true/false 

Write  time-averaged  data  if 
true 

Write  snapshots  of  data  if 
false 

=  .false. 

restart 

true/false 

Initialize  using  restart  file 

=  .true. 

print_points 

true/false 

Print  diagnostic  data  for 
two  grid  points 

=  .true. 

print  global 

true/false 

Print  diagnostic  data,  global 
sums 

=  .true. 

kitd 

0/1 

If  0,  delta  function  ITD 
approx. 

If  1 ,  linear  remapping  ITD 
approx. 

=  1 

kcatbound 

0/1 

If  0,  original  category 
boundary  formula 

If  1 ,  new  category  boundary 
formula 

=  1 

kdyn 

0/1 

If  0,  EVP  dynamics  OFF 

If  1,  EVP  dynamics  ON 

=  1 

kstrength 

0/1 

If  0,  [15]  ice  strength 
formulation 

If  1,  [35]  ice  strength 
formulation 

=  1 

krdg_partic 

0/1 

If  0,  old  ridging 
participation  function 

If  1 ,  new  ridging 
participation  function 

=  1 

krdgredist 

0/1 

If  0,  old  ridging 
redistribution  function 

If  1 ,  new  ridging 
redistribution  function 

=  0 

evpdamping 

true/false 

If  true,  damp  elastic  waves, 
[17] 

=  .false. 

advection 

remap,  mpdata, 
upwind 

remap:  Linear  remapping 

advection 

mpdata:  2nd  order 

=  ’remap’ 
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Name 

Type/Options 

Description 

Default  Values  /  Directory 
Location 

MPDATA 

upwind:  1st  order  MPDATA 

grid_type 

rectangular, 

displaced_pole 

rectangular:  Defined  in 
rectgrid 

displaced_pole:  Read  from 
file  in  popgrid 

=  ’pips’ 

gridfile 

filename 

Name  of  grid  file  to  be  read 

=  ’grid_cice_1280x720.r’ 

kmt  file 

filename 

Name  of  land  mask  file  to 
be  read 

=  'kmt' 

dumpfile 

filename  prefix 

Output  file  for  restart  dump 

=  ’iced’ 

restart  dir 

path/ 

path  to  restart  directory 

=  ’/scr/posey/pips3c/’ 

pointer  file 

pointer  filename 

Contains  restart  filename 

’/scr/posey/pips3c/ice. restart  file’ 

hist  dir 

path/ 

path  to  history  output 
directory 

=  ’/scr/posey/pips3c/’ 

history  file 

filename  prefix 

Output  file  for  history 

=  ’iceh’ 

diagfile 

filename 

Diagnostic  output  file 

=  ’ice_diag.d’ 

oceanmixed  ice 

true/false 

Active  ocean  mixed  layer 
calculation 

=  .true. 

albicev 

0  <  a  <  1 

Visible  ice  albedo  for 
thicker  ice 

=  0.65 

albicei 

0  <  a  <  1 

Near  infrared  ice  albedo  for 
thicker  ice 

=  0.65 

albsnowv 

0  <  a  <  1 

Visible,  cold  snow  albedo 

=  0.85 

albsnowi 

0  <  a  <  1 

Near  infrared,  cold  snow 
albedo 

=  0.85 

ycycle 

integer 

No.  of  years  in  forcing  data 
cycle 

=  1 

fyear  init 

yyyy 

First  year  of  atmospheric 
forcing  data 

=  2008 

atm  data  dir 

path/ 

Path  to  atm  forcing  data 
directory 

=  ’/scr/posey/pips3c/data_in/’ 

ocn  data  dir 

path/ 

Path  to  oceanic  forcing  data 
directory 

=  7scr/posey/pips3c/data_in/’ 

Table  1:  Ice  namelist  parameters. 
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4.1. 3.2  Ice  Fields  Namelist  (icefields  nml) 


Name 

Description 

Default  Values 

f  hi 

Ice  thickness 

=  .true. 

f  hs 

Snow  thickness 

=  .true. 

fTsfc 

Temperature  of  ice/snow 
top  surface  (in  category  n) 

=  .false. 

f  aice 

Ice  concentration 

=  .true. 

f_uvel 

x-componcnt  of  velocity 

=  .true. 

f_vvel 

v-component  of  velocity 

=  .true. 

f_fswdn 

Incoming  shortwave 
radiation  down 

=  .false. 

fflwdn 

Incoming  longwave 
radiation  down 

=  .fasle. 

f_snow 

Snowfall  rate 

=  .false. 

f  snow  ai 

Snowfall  rate  weighted  by 
aice 

=  .false. 

f  rain 

Rainfall  rate 

=  .false. 

f  rain  ai 

Rainfall  rate  weighted  by 
aice 

=  .false. 

f_sst 

Sea  surface  temperature 

=  .true. 

f_sss 

Sea  surface  salinity 

=  .true. 

fuocn 

Ocean  current,  x  direction 

=  .true. 

fvocn 

Ocean  current,  y  direction 

=  .true. 

f  frzmlt 

Freezing/melting  potential 

=  .false. 

f_fswabs 

Absorbed  shortwave 
radiation 

=  .false. 

f  fswabs  ai 

Absorbed  shortwave 
radiation  weighted  by  aice 

=  .false. 

f  albsni 

Snow/ice  broad  band 
albedo 

=  .false. 
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Name 

Description 

Default  Values 

falvdr 

=  .false. 

f  alidr 

=  .false. 

f  flat 

Latent  heat  flux 

=  .false. 

f  flat  ai 

Latent  heat  flux  weighted 
by  aice 

=  .false. 

f  fsens 

Sensible  heat  flux 

=  .false. 

f  fsens  ai 

Sensible  heat  flux  weighted 
by  aice 

=  .false. 

fflwup 

Incoming  longwave 
radiation  upward 

=  .false. 

fflwupai 

Incoming  longwave 
radiation  upward  weighted 
by  aice 

=  .false. 

f_evap 

Evaporation  water  flux 

=  .false. 

fevapai 

Evaporation  water  flux 
weighted  by  aice 

=  .false. 

fTref 

2m  atmospheric  reference 
temperature 

=  .false. 

f  Qref 

2  m  atmospheric  reference 
specific  humidity 

=  .false. 

fcongel 

Basal  ice  growth 

=  .false. 

f  frazil 

Frazil  ice  growth 

=  .false. 

fsnoice 

Snow-ice  fonnation 

=  .false. 

f  meltt 

Top  ice  melt 

=  .false. 

f  meltb 

Basal  ice  melt 

=  .false. 

f  meltl 

Lateral  ice  melt 

=  .false. 

f  fresh 

Fresh  water  flux  to  ocean 

=  .false. 

f  fresh  ai 

Fresh  water  flux  to  ocean 
weighted  by  aice 

=  .false. 

f  fsalt 

Net  salt  flux  to  ocean 

=  .false. 

f  fsalt  ai 

Net  salt  flux  to  ocean 
weighted  by  aice 

=  .false. 
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Name 

Description 

Default  Values 

f  fhnet 

Net  heat  flux  to  ocean 

=  .false. 

f  fhnet  ai 

Net  heat  flux  to  ocean 
weighted  by  aice 

=  .true. 

f_fswthru 

Shortwave  penetrating  to 
ocean 

=  .false. 

f  fswthru  ai 

Shortwave  penetration  to 
ocean  weighted  by  aice 

=  .false. 

f  strairx 

Stress  on  ice  by  air  in  the 
x-direction  (centered  in  U 
cell) 

=  .true. 

f  strairy 

Stress  on  ice  by  air  in  y- 
direction  (centered  in  T 
cell) 

=  .true. 

f  strtltx 

Surface  stress  due  to  sea 
surface  slope  in  x-direction 

=  .false. 

f  strtlty 

Surface  stress  due  to  sea 
surface  slope  in  y-direction 

=  .false. 

f  strcorx 

Coriolis  stress  (x) 

=  .false. 

f  strcory 

Coriolis  stress  (y) 

=  .false. 

f  strocnx 

Ice-ocean  stress,  x  dir.  (U 
cell) 

=  .true. 

f  strocny 

Ice-ocean  stress,  y-dir.  (T 
cell) 

=  .true. 

f  strintx 

Divergence  of  internal  ice 
stress,  x-direction 

=  .true. 

f  strinty 

Divergence  of  internal  ice 
stress,  y-direction 

=  .true. 

f  strength 

Ice  strength  (pressure) 

=  .true. 

fopening 

Lead  area  opening  rate 

=  .true. 

fdivu 

Strain  rate  I  component, 
velocity  divergence 

=  .false. 

f  shear 

Strain  rate  II  component 

=  .false. 

f_sigl 

Principal  stress 
components  (diagnostic) 

=  .false. 
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Name 

Description 

Default  Values 

f_sig2 

Principal  stress 
components  (diagnostic) 

=  .false. 

fdvidtt 

Ice  volume  tendency  due  to 
thennodynamics 

=  .false. 

fdvidtd 

Ice  volume  tendency  due  to 
dynamic  s/ transport 

=  .false. 

fdaidtt 

Ice  area  tendency  due  to 
thennodynamics 

=  .false. 

fdaidtd 

Ice  area  tendency  due  to 
dynamic  s/ transport 

=  .false. 

f  mlt  onset 

Day  of  year  that  surface 
melt  begins 

=  .false. 

f  frz  onset 

Day  of  year  that  freezing 
begins 

=  .false. 

fdardgldt 

Ice  area  ridging  rate 

=  .false. 

f_dardg2dt 

Ridge  area  formation  rate 

=  .false. 

fdvirdgdt 

Ice  volume  ridging  rate 

=  .false. 

f  hisnap 

Ice  volume  snapshot 

=  .false. 

f  aisnap 

Ice  area  snapshot 

=  .false. 

f  aicel  (through 

5) 

Ice  concentration  in  grid 
cell  in  categories  1  through 

5 

=  .true. 

f  aice6  (through 
10) 

Ice  concentration  in  grid 
cell  in  categories  6  through 
10 

=  .false. 

f  vicel  (through 

5) 

Volume  per  unit  area  of  ice 
in  categories  1  through  5 

=  .true. 

f_vice6 
(through  10) 

Volume  per  unit  area  of  ice 
in  categories  6  through  10 

=  .false. 

Table  2:  Ice  field  namelist  parameters. 


4. 1. 4  PIPS  3. 0  Macros  File 

Macros.AIX.txt  is  a  file  containing  macros  necessary  for  compiling  the  PIPS  3.0  code  on  the 
NAVOCEANO  IBM  platfonn  “Babbage”.  A  detailed  description  of  the  macro  is  provided  in  the 
PIPS  3.0  User’s  Manual  [2].  The  primary  macros  utilized  in  this  file  are  summarized  below. 
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Macro 

Description 

-lmass 

IBM  -  tuned  intrinsic  library 

-qsmp=noauto 

enables  SMP  directives,  but  does  not  add  any 

-qstrict 

don't  turn  divides  into  multiplies,  etc 

-qhot 

higher-order  -transformations  (eg.  loop  padding) 

-qalias=noaryoverlp 

assume  no  array  overlap  with  respect  to  equivalence,  etc 

-qmaxmem=l 

memory  available  to  compiler  during  optimization 

-qipa=level=2 

InterProcedure  Analysis  (eg.  inlining)  =>  slow  compiles 

l 

P> 

i 

P 

enable  profiling  (use  in  both  FFLAGS  and  LDFLAGS) 

-qreport 

for  smp/omp  only 

-bmaxdata : 0x80000000 

use  maximum  allowed  data  segment  size 

-g 

always  leave  it  on  because  overhead  is  minimal 

-qf  lttrap=... 

enable  default  sigtrap  (core  dump) 

-C 

runtime  array  bounds  checking  (runs  slow) 

-qini  tauto=... 

initializes  automatic  variables 

Table  3:  PIPS  3.0  macros  from  Macros.AIX.txt  file. 


4. 1. 5  PIPS  3. 0  Makefile 

A  makefile  is  used  in  the  execution  of  PIPS  3.0  at  NAVOCEANO.  A  comprehensive  description 
of  the  makefile  is  provided  in  the  PIPS  3.0  User’s  Manual  [2].  The  command  line  variables  and 
usage  examples  are  defined  below. 

Command-line  variables: 

1 .  MACFILE=<file>  ~  The  macros  definition  file  to  use/include  in  a  run. 

2.  EXEC=<name>  ~  The  name  given  to  an  executable.  The  default  is  a. out. 

3.  VPATH=<vpath>  ~  VPATH,  default  is  .  (cwd  only). 

4.  SRCS=<files>  ~  A  list  of  source  files.  The  default  is  all  .c  .F  .F90  files  in  VPATH. 

5.  VPFILE=<file>  ~  A  file  with  a  list  of  directories.  It  is  used  to  create  VPATH. 

6.  SRCFILE=<file>  ~  A  file  with  a  list  of  source  files.  It  is  used  to  create  SRCS. 

7.  DEPGEN=<exec>  ~  A  dependency  generator  utility,  with  a  default  of  makdep. 

8.  <macro  defns>  ~  Any  macro  definitions  found  in  this  file  or  the  included  MACFILE  will 
be  over-ridden  by  command  line  macro  definitions. 

9.  MODEL=<model>  ~  A  standard  macro  definition,  often  found  in  the  included 
MACFILE.  It  is  used  to  trigger  special  compilation  flags. 

Usage  examples: 

%  gmake  MACFILE=Macros.AIX  VPFILE=Filepath  MODEL=ccm3  EXEC=atm 
%  gmake  MACFILE=Macros.AIX  VPFILE=Filepath  SRCFILE=Srclist  EXEC=pop 

14 


PIPS  3.0  SDD 


%  gmake  MACFILE=Macros.C90  VPATH="dirl  dir2M  SRCS="filel.c  file2.F90" 
%  gmake  MACFILE=Macros.SUN  SRCS="test.F" 

4.2  PIPS  3.0  Software  Organization  and  Implementation 

4.2.1  Logical  Component  Call  Trees 

4. 2. 1.1  Primary  Tree  (beginning  at  the  program  ‘icemodel ) 

icemodel 

I 

+-ice_albedo — t — ice_kinds_mod 
I  I 

I  H — ice_domain — I — ice_kinds_mod 

I  I 

I  H — ice  model  slze--ice  kinds  mod 

I 

+-ice_calendar--ice_constants — I — ice_kinds_mod 
I  '  '  I 

I  H — ice_domain 

I 

+-ice_coupling — I — ice_kinds_mod 

I  I 

I  h — ice_model_size 

I  { 

I  +-ice_constants 

I  I 

I  +-ice_calendar 

I  I 

I  +-ice_state-+-ice_kinds_mod 

I  I  I 

I  I  h — ice_model_size 

I  I  I 

I  I  H — ice_domain 

I  I 

I  h — ice_flux — l — ice_kinds_mod 

I  I  I 

I  I  H — ice_domain 

I  I  I 

I  I  +-ice_constants 

I  I 

I  +-ice_albedo 

I  I 

I  H — ice_mpi_internal — I — ice_kinds_mod 

I  I  ~  _  I 

I  I  H — ice_domain 

I  I 

I  H — ice_timers — I — ice_kinds_mod 

I  I  ""  I 

I  I  +-ice_constants 

I  I 

I  h — ice_f ileunits--ice_kinds_mod 

I  I 

I  H — ice_work — I — ice_kinds_mod 

I  _  I 

I  H — ice_domain 

I 

+-ice_diagnostics — I — ice^domain 
I  r  I 

I  +-ice_constants 

I  I 

I  +-ice  calendar 
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I  h — ice_f ileunits 

I  I 

I  H — ice_work 

I 

H — ice_domain 
I 

H — ice_dyn_evp — I — ice_kinds_mod 
I  “  I 

I  H — ice_domain 

I  I 

I  +-ice_constants 

I  I 

I  +-ice_state 

I  I 

I  H — ice_work 

I 

+-ice_f ileunits 
I 

H — ice_f lux_in — i--ice_kinds_mod 
I  ~~  _  I 

I  H — ice  domairi 

I  I 

I  +-ice_constants 

I  I 

I  H — ice_flux 

I  I 

I  +-ice_calendar 

I  I 

I  h — ice_read_write — t--ice_model_size 

I  I  _  “  I 

I  I  H — ice_domain 

I  I  I 

I  I  H — ice_mpi_internal 

I  I  t 

I  I  +-ice_f ileunits 

I  I  I 

I  I  H — ice_work 

I  I 

I  H — ice_f ileunits 

I 

H — ice_grid — i--ice_kinds_mod 

C  '  "  '  ’  :f 

I  +-ice_constants 

I  I 

I  H — ice_domain 

I  I 

I  H — ice_f ileunits 

I  I 

I  H — ice_mpi_internal 

I  I 

I  H — ice_work 

I 

H — ice_hi story — I — ice_kinds_mod 

f  I 

I  H — ice_domain 

I  I 

I  +-ice_read_write 

I  I 

I  H — ice_f ileunits 

I  I 

I  H — ice_work 

I 

H — ice_init--ice_domain 
I 

H — ice_itd-H — ice_kinds_mod 
I  _  I 

|  H — ice  model  size 
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I  +-ice_constants 

I  I 

I  +-ice_state 

I  I 

I  H — ice_f ileunits 

I 

h — ice_itd_l inear — I — ice_model_size 
I  I 

I  H — ice_kinds_mod 

I  I 

I  H — ice_  domairi 

I  I 

I  +-ice_constants 

I  I 

I  +-ice_state 

I  I 

I  +-ice_itd 

I  I 

I  +-ice_calendar 

I  I 

I  H — ice_f ileunits 

I 

H — ice_kinds_mod 
I 

h — ice_mechred — I — ice_model_size 

T-  ! 

I  +-ice_constants 

I  I 

I  +-ice_state 

I  I 

I  +-ice_itd 

I  I 

I  +-ice_grid 

I  I 

I  H — ice_f ileunits 

I  I 

I  H — ice_domain 

I  I 

I  +-ice  calendar 


I  H — ice_work 

I 

H — ice_mpi_internal 
I 

H — ice_ocean — I — ice_kinds_mod 
I  I 

I  +-ice_constants 

I 

+  -ice_scaling — I — ice_domain 

(  I 

I  H — ice_kinds_mod 

I  I 

I  +-ice_constants 

I  I 

I  +-ice_state 

I  I 

|  +-ice_flux 

I  I 

I  +-ice_grid 

I 

H — ice_therm_vertical — I — ice_model_size 
I  I 

I  H — ice_kinds_mod 

I  I 

I  H — j. ce  domairi 

I  I 

I  H — ice  fileunits 


I  +-ice_constants 

I  I 

I  +-ice_calendar 

I  I 

I  +-ice_grid 

I  I 

I  +-ice_state 

I  I 

I  H — ice_flux 

I  I 

I  +-ice_itd 

I  I 

I  +-ice_diagnostics 

I 

H — ice_therm_itd — I — ice_kinds_mod 
I  ~~  ~~  I 

I  h — ice_model_size 

I  f 

I  +-ice_constants 

I  I 

I  H — ice_domain 

I  I 

I  +-ice_state 

I  I 

I  H — ice_flux 

I  I 

I  +-ice_diagnostics 

I  I 

I  +-ice_calendar 

I  I 

I  +-ice_grid 

I  I 

I  +-ice_itd 

I 

+-ice_timers 

I 

h — ice_transport_mpdata — I — ice_model_size 
'[:  ^  I 

I  H — ice_domain 

I  I 

I  +-ice_constants 

I  I 

I  +-ice_grid 

I  I 

I  h — ice_f ileunits 

I  I 

I  +-ice_init 

I  I 

I  H — ice_work 

I 

h — ice_transport_remap — I — ice_model_size 
f  ~  I 

I  H — ice_kinds_mod 

I  I 

I  H — ice_domain 

I  I 

I  +-ice_constants 

I  I 

I  +-ice_grid 

I  I 

I  H — ice_f ileunits 

I  I 

I  +-ice_calendar 

I  I 

I  +-ice_state 

I  I 

I  +-ice_timers 

I  I 

I  +-ice  itd 
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I  H — ice_work 

I 

+  - (shr_msg_stdio) 

I 

H — setup_mpi — l--ice_mpi_internal 
I  ] 

I  +-ice_coupling 

I  ] 

I  +-ice_coupling_setup-+- (cpl_interface_init) 

II  I 

I  I  +- (shr_sys_flush) 

I  I 

|  H —  (mpi_init) 

I  I 

|  H —  (mpi_comm_dup) 

I  I 

|  h —  (mpi_comm_size) 

I  I 

I  H —  (mpi_comm_rank) 

I  I 

I  +- (abort_ice) 

I  I 

|  +- (mpi_type_vector ) 

I  I 

|  +- (mpi_type_commit) 

I 

h — ice_timer_clear 

I 

H — ice_timer_s tart- -timers — I —  (mpi_wtime) 

I  I 

I  +- (irtc_rate) 

I  I 

I  H —  (rtc) 

I 

+-init_constants 

I 

+-input_data-+-ice_albedo 
I  I 

I  +-ice_diagnostics 

I  I 

I  +-ice_history 

I  I 

I  +-ice_calendar 

I  I 

I  +-ice_dyn_evp 

I  i 

I  H — ice_flux_in 

I  I 

I  +-ice_coupling 

I  I 

I  +-ice_bcast_iscalar-- (mpi_bcast) 

I  I 

I  +- (abort_ice) 

I  I 

I  +-ice_bcast_rscalar-- (mpi_bcast) 

I  I 

I  +-ice_bcast_char-- (mpi_bcast) 

I  I 

I  +-ice_bcast_logical-- (mpi_bcast) 

I 

H — init_grid 
I  I 

I  +-global_scatter-+-ice_model_size 

II  I 

I  |  +-ice_constants 

II  I 

|  |  +- (mpi_irecv) 
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|  |  +- (mpi_isend) (mpi_sendrecv) 

II  I 

|  |  +- (mpi_irecv) 

II  I 

|  |  +- (mpi_isend) 

II  I 

|  |  +- (mpi_wait) 

I  I 

I  +-popgrid-+-ice_read_write 

II  I 

I  |  +-ice_open 

II  I 

I  |  +-ice_read-+-global_scatter 

II  I 

I  |  +-ice_bcast_logical 

I  I 

I  H — pipsgrid — I — ice_read_write 

II  I 

I  |  +-ice_open 

II  I 

I  |  +-ice_read 

I  I 

I  +-columngrid-+-global_scatter 

II  I 

I  |  H — ice_model_size 

II  I 

I  |  +- (abort_ice) 

I  I 

I  +-rectgrid-+-global_scatter 

II  I 

I  |  H — ice_model_size 

I  I 

I  +-bound — bound_ijn — I — ice_timers 

II  ~  I 

I  |  +-ice_timer_start 

I  I  I 

|  |  +-+-ice_timer_stop--timers 

I  I 

I  +- (abort_ice) 

I  I 

I  +-tlatlon — i — global_gather — I — ice_model_size 

II  I  ~~  I 

I  |  I  +-ice_constants 

II  I  I 

||  |  +- (mpi_irecv) 

II  I  1 

||  |  +- (mpi_barrier ) 

II  I  I 

||  |  +- (mpi_isend) 

II  I  I 

||  |  +- (mpi_wait) 

II  I  I 

|  |  |  +- (mpi_waitall) 

III 

I  |  +-global_scatter 

II  I 

I  |  h — ice_model_size 

II  I 

I  |  +-bound 

I  I 

|  H — makemask--bound 
I 

H — init_remap — I — bound 
I  I 

I  +- (abort_ice) 

I 

H — init  calendar 
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H — init_hist — i — ice_bcast_i scalar 

i  I 

I  +-ice_bcast_logical 

I  I 

I  +-ice_constants 

I  I 

I  +-ice_calendar 

I  I 

I  H — ice_flux 

I  I 

I  +- (abort_ice) 

I  I 

I  +- (shr_sys_flush) 

I 

h — init_evp — I — ice_calendar 
I  _  I 

I  h — ice_f ileunits 

I  I 

I  H — ice_flux 

I  I 

I  +- (ulat) 

I  I 

I  +- (bound) 

I 

H — init_f lux — i — ice_constants 
I  ] 

I  H — ice_flux 

I  I 

|  +-init_f lux_atm--ice_state 

I  1 

|  H — init_f lux_ocn 

I 

H — init_thermo_vertical--ice_itd 
I 

+-init_mechred 

I 

+-init_itd-- (abort_ice) 

I 

h — calendar- -ice_f ileunits 
I 

H — init_cpl 
I  I 

I  +-calendar 
I  I 

I  +-ice_bcast_iscalar 
I  I 

I  +-ice_bcast_rscalar 

I  I 

I  +- (shr_sys_flush) 

I  I 

I  +-  (tlon) 

I  I 

I  +-(tlat) 

I  I 

I  +- (tarea) 

I  I 

|  +- (hm) 

I  I 

I  +- (cpl_interface_contractinit) 

I  I 

I  +- (cpl_interface_ibufrecv) 

I  I 

I  +-to_coupler 

I  I 

I  H — get_sum--ice_global_real_sum-- (mpi_allreduce) 

I  I 

I  +-ice  timer  start 
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I  +-ice_timer_stop 

II 

I  +- (anglet) 

I  I 

I  +- (tmask) 

I  I 

I  +- (shr_sys_flush) 

II 

I  +- (cpl_interface_contractsend) 

I  I 

I  +-  (tarea) 

I  I 

I  +- (bound) 

I 

h — init_getf lux — I — ncar_f iles--f ile_year 

f  I 

I  H — sss_clim — l — ice_work 

I  (  “  j 

I  I  +-ice_open 

I  I  I 

I  I  +-ice_read 

I  I  I 

|  |  +-complete_getf lux_ocn 

I  I 

I  +-sst_ic-+-ice_open 

I  I 

I  +-ice_read 

I 

h — init_state — I — ice_model_size 

I  I 

I  +-ice_constants 

I  I 

I  H — ice__flux 

I  I 

I  +-ice_grid 

I  I 

I  +-ice_state 

I  I 

I  +-bound 

I  I 

|  H —  (hin_max) 

I  I 

I  +- (ilyrl) 

I  I 

I  +- (tmlt) 

I  I 

I  +- (aggregate) 

I  I 

I  +- (bound_aggregate) 

I 

H — restartf  ile — i--ice_model_size 

I  I 

I  H — ice_mpi_internal 

I  I 

I  +-ice_open 

I  I 

I  +-ice_bcast_iscalar 

I  j. 

I  +-ice_bcast_rscalar 

I  I 

I  +-ice_read 

I  I 

I  H — ice_flux 

I  I 

I  +-ice_grid 

I  I 

I  +-ice_calendar 

I  I 

I  +-ice  state 
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I  +-ice_dyn_evp 

I  I 

I  +-ice_coupling 

I  I 

I  +-lenstr 

I  I 

I  +- (bound_state) 

I  I 

I  +-bound 

I  I 

I  +- (aggregate) 

I  I 

I  +- (bound_aggregate) 

I 

+-albedos-+-ice_constants 
I  I 

I  +-ice_state 

I  I 

I  +- (tmask) 

I 

H — init_diags — I — ice_mpi_internal 
I  \ 

I  +-global_gather 

I  I 

I  +-(tlat_g) 

I  I 

I  +-(tlon_g) 

I 

+-init_diagnostics--ice_state 

I 

h — f rom_coupler — I — get_sum 
I  “  I 

I  +-ice_timer_start 

I  I 

I  +-ice_timer_stop 

I  I 

I  +- (cpl_interface_contractrecv) 

I  I 

I  +- (tarea) 

I  I 

|  +-  (hm) 

I  I 

I  +- (bound) 

I  I 

I  +- (anglet) 

I  I 

I  +- (t2ugrid) 

I 

+  -getf  lux-H — sss_sst_restore — I — interp_coef f_monthly 

I  I  ’  I 

I  I  +-read_clim_data-+-ice_open 

II  I  I 

I  I  I  +-ice_read 

II  I  I 

I  I  I  +-ice_diagnostics 

II  I 

I  I  +-interpolate_data 

I  I  I 

|  |  +-complete_getf lux_ocn 

II 

I  H — ncar_bulk_dat — I — interp_coef f_monthly 

I  I  ""  _  I 

I  I  +-read_clim_data 

II  I 

I  I  +-interpolate_data 

II  I 

I  I  +-interp_coef f 
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I  I  +-read_data-+-ice_open 

II  [ 

I  I  +-ice_read 

II  I 

I  I  +-ice_diagnostics 

II  r 

|  |  +-file_year 

II 

I  +-prepare_forcing--ice_state 

I 

+  -init_mass_diags — I — ice_mpi_internal 

t  I 

I  +-ice_state 

I  I 

I  +-get_sum 

I  I 

I  +-global_gather 

I 

H — mixed_layer — i — ice_f lux 

’  '  I 

I  +-ice_calendar 

I  I 

I  +-ice_grid 

I  I 

I  +-ice_state 

I  I 

I  +-ice_albedo 

I  I 

I  +- (atmo_boundary_layer ) 

I 

H — thermo_vertical 
I  I 

I  H — ice_work 

I  I 

|  H — init_f lux_atm 

II 

|  +-init_diagnostics 
I  I 

|  +-merge_f luxes--ice_state 

I  I 

I  +-ice_timers 

I  I 

I  +-ice_ocean 

I  I 

I  +-ice_timer_start 

I  I 

|  H — f r zmlt_bottom_lateral 

II 

|  +-init_thermo_vars 

I  I 

I  +- (atmo_boundary_layer ) 

I  I 

|  H — init_vertical_prof ile — I — print_state — I — ice_model_size 

|  I  ^  I  ^  1 

I  |  I  H — ice_kinds_mod 

II  I  ( 

I  |  I  +-ice_state 

II  I  I 

I  |  I  +-ice_flux 

II  I  I 

I  |  I  +- (ilyrl) 

I  I  I 

I  |  +- (abort_ice) 

I  I 

I  +-temperatur exchanges -+-print_state 

II  I 

I  |  +-conductivity 
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I  |  +-absorbed_solar--ice_albedo 

I  I  I 

I  |  +-surface_f luxes 

II  I 

I  |  +-tridiag_solver 

II  I 

I  |  +- (abort_ice) 

I  I 

|  H — thickness_changes 
I  I 

I  +-conservation_check_vthermo — i — print_state 

II  ~~  I 

I  |  +- (abort_ice) 

I  I 

I  H — add_new_snow 

I  I 

I  +-update_state_vthermo 

II 

I  +-ice_timer_stop 
I 

+-scale_f luxes--ice_albedo 

I 

+-to_coupler 

I 

H — thermo_itd — I — aggregate — I — ice_domain 
I  t  I 

I  I  +-ice_flux 

I  I  I 

I  I  +-ice_grid 

I  I 

|  H — init_f lux_ocn 

I  I 

I  +-reduce_area--ice_grid 

I  t 

I  H — rebin  (44) — I — ice_grid 

I  I  I 

|  |  H — shift_ice — I — ice_flux 

I  I  _  I 

I  I  H — ice_work 

II  I 

I  I  +- (abort_ice) 

I  I 

I  H — zap_small_areas — I — ice_flux 

I  I  ""  _  I 

I  I  +-ice_calendar 

I  I  I 

I  I  +- (abort_ice) 

I  I 

I  +-ice_timers 

I  I 

I  H — ice_itd_linear 

I  I 

I  H — ice_therm_vertical 

I  I 

I  +-ice_timer_start 

I  I 

I  +-ice_timer_stop 

I  I 

|  +-linear_itd-+-aggregate_area-- (abort_ice) 

I  I  I 

I  I  H — column_sum 

I  I  I 

I  I  H — column_conservation_check-- (abort_ice) 

I  I  I 

|  |  +-shift_ice 

I  I  I 

|  |  +-fit_line 

I  I 

I  H — add  new  ice — I — column  sum 
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I  I  H — column_conservation_check 

II 

I  +-lateral_melt 

I  I 

I  +-freeboard 

I 

H — evp — l — ice_timers 
I  I 

I  +-ice_timer_start 

I  I 

I  H — evp_prep — I — ice_flux 

I  [  '  '  1 

I  I  +-ice_calendar 

I  I  I 

I  I  +- (bound) 

I  I  I 

I  I  +- (tmask) 

I  I  I 

I  I  +- (t2ugrid) 

I  I  I 

I  I  +- (to_ugrid) 

II  I 

I  I  +- (iceumask) 

I  I  I 

I  I  +- (ice_strength) 

II  t. 

I  I  +-(icetmask) 

I  I  I 

I  I  +- (umask) 

I  I 

I  +-stress-+- (cyp) 

I  I  I 

I  I  H —  (dyt) 

I  I  I 

I  I  +-(cxp) 

I  I  I 

I  I  +-(dxt) 

I  I  I 

I  I  +- (cym) 

I  I  I 

I  I  +- (cxm) 

I  I  I 

I  I  +- (tarear) 

I  I  I 

I  I  +- (tinyarea) 

I  I 

I  +-stepu-+-ice_f lux 

II  I 

I  I  H — (dxt2 ) 

I  I  I 

I  I  +-(dyt2) 

I  I  I 

I  I  +-(dyt4) 

I  I  I 

I  I  +- (dxhy) 

I  I  I 

I  I  +- (dyhx) 

I  I  I 

I  |  +-(dxt4) 

I  I  I 

I  I  +- (bound_narr_ne) 

II  I 

I  I  +- (iceumask) 

I  I  I 

I  I  +- (uarear) 

I  I 

I  +- (bound_sw) 
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I  H — evp_finish — I — ice_flux 

I  j  I 

I  I  +- (iceumask) 

I  I  I 

I  I  +- (u2tgrid) 

I  I 

I  +-ice_timer_stop 

I 

H — transport_remap — I — bound 
I  ~~  I 

I  H — bound_sw--bound_i j n 

I  I 

I  +-bound_state-+-ice_grid 

I  I  I 

I  I  H — bound_narr--bound_i j n 

I  I 

I  H — bound_narr 

I  I 

I  +-ice__timer_start 

I  I 

I  +-ice_timer_stop 

I  I 

I  +-departure_points 

I  I 

I  +-locate_triangles 

I  I 

I  +-triangle_coordinates 

I  I 

|  H — make_masks 

I  I 

I  H — conserved_sums — I — ice_mpi_internal 

I  I  I 

I  I  +-ice_global_real_sum 

I  I 

I  H — construct_f ields--limited_gradient 

I  I 

|  H — f lux_integrals 

I  I 

I  +-update_f ields-- (abort_ice) 

I  I 

I  +-global_conservation-- (abort_ice) 

I  I 

I  +-load_tracers 

I  I 

I  H — local_max_min--bound 

I  I 

I  +-check_monotonicity-- (abort_ice) 

I  I 

|  +-unload_tracers 

I 

H — transport_mpdata — I — bound_narr 
I 

I  +-ice_flux 

I  I 

I  +-ice_timers 

I  I 

I  +-ice_state 

I  I 

I  +-ice_itd 

I  I 

I  +-ice_timer_start 

I  I 

I  +-check_state--ice_f lux 

I  I 

|  +-mpdata-+-bound 

I  I  I 

I  I  +-bound  narr 
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I  I  +-ice_calendar 

I  I  I 

I  I  +-ice_state 

I  ii 

I  I  +- (abort_ice) 

I  I 

I  +-ice_timer_stop 

I 

h — ridge_ice — i — ice_timers 
I  _  I 

I  H — ice_flux 

I  I 

I  +-ice_timer_start 

I  I 

I  H — ridge_prep--asum_ridging 

I  I 

I  h — ridge_shift — I — column_sum 

I  I  _  I 

I  I  +-column_conservation_check 

II  I 

I  I  +- (abort_ice) 

I  I 

I  H — asum_ridging 

I  I 

I  +- (abort_ice) 

I  I 

I  +-ice_timer_stop 

I 

h —  zap_small_areas 

I 

h — rebin 

I 

+-aggregate 

I 

+-scale_hist_f luxes 

I 

h — runtime_diags — I — ice_model_size 
I  _  I 

I  H — ice_flux 

I  I 

I  +-ice_albedo 

I  I 

I  H — ice_mpi_internal 

I  I 

I  +-ice_state 

I  I 

I  +- (mask_n) 

I  I 

I  +- (mask_s) 

I  I 

I  I — ice_global_real_maxval-- (mpi_allreduce) 

I  I 

I  +-get_sum 

I  I 

I  +- (ulat) 

I  I 

I  +- (ulon) 

I  I 

I  +- (opening) 

I  I 

I  +- (dardgldt) 

I  I 

I  +- (strintx) 

I  I 

I  +- (strinty) 

I  I 

I  +- (athorn) 
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I  H — ice_global_real_sum 

I  I 

I  +-global_gather 

I  I 

I  +- (shr_sys_flush) 

I 

h — ice_write_hist — I — ice_f lux 
I  ~~  ~~  I 

I  +-ice_albedo 

I  I 

I  +-ice_grid 

I  I 

I  +-ice_calendar 

I  I 

I  +-ice_state 

I  I 

I  +-ice_dyn_evp 

I  I 

I  +-ice_constants 

I  I 

I  +- (opening) 

I  I 

I  +- (dardgldt) 

I  I 

I  +- (dardg2dt) 

I  I 

I  +- (dvirdgdt) 

I  I 

|  +-principal_stress 

I  I 

I  +-icecdf-+-global_gather 

I  I 

I  H — ice_model_size 

I  I 

I  h — ice_mpi_internal 

I  I 

I  +-ice_constants 

I  I 

I  +-ice_grid 

I  I 

I  +-ice_calendar 

I  I 

I  +-lenstr 

I  I 

I  +- (nf_create) 

I  I 

|  H —  (nf_def_dim) 

I  I 

I  +- (nf_def_var) 

I  I 

I  +- (nf_put_att_text) 

I  I 

I  +- (nf_put_att_real) 

I  I 

I  +- (nf_enddef ) 

I  I 

|  H —  (nf_inq_varid) 

I  I 

I  +- (nf_put_var_real) 

I  I 

I  +- (nf_put_vara_real) 

I  I 

I  +- (nf_close) 

I 

H — dumpf ile — I — ice_model_size 
I  I 

I  +-ice_open 
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I  +-ice_write--global_gather 

I  I 

I  H — ice_flux 

I  I 

I  +-ice_grid 

I  I 

I  +-ice_calendar 

I  I 

I  +-ice_state 

I  I 

I  +-ice_dyn_evp 

I  I 

I  +-ice_coupling 

I  I 

I  +-lenstr 

I 

+-ice_timer_stop 

I 

H —  ice_time  r _p  rint  — I —  ice_domain 
I  I 

I  H — ice_mpi_internal 

I  I 

I  H — ice_f ileunits 

I  I 

I  h — ice_global_real_minval-- (mpi_allreduce) 

I  I 

I  h — ice_global_real_maxval 

I 

H — end_run — (mpi_f inalize) 

I 

+-exit_coupler-+- (mpi_abort) 

I 

H —  (cpl_interface_f inalize) 

4.2. 1.2  Detached  Call  Trees  (for  extraneous  modules  and  subroutines) 

1.  ice_exit--ice_kinds_mod 

2.  atmo_boundary_layer--ice_grid 

3 .  ice_global_real_minmax — I — ice_f ileunits 

I 

i — ice_global_real_minval 

I 

H — ice_global_real_maxval 

4.  bound_aggregate-+-ice_grid 

I 

+-bound 

5.  bound_narr_ne--bound_i j n 

6.  u2tgrid-+-bound 

I 

+-to_tgrid 

7.  ice_atmo — I — ice_domain 

I 

+-ice_constants 

I 

H — ice_f lux 
I 

+-ice  state 
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8.  debug_ice — I — ice_kinds_mod 

I 

+-ice_itd 

I 

+-ice_diagnostics 

I 

+-ice_mpi_internal 

I 

+-ice_grid 

I 

+-print_state 

9.  ice_strength--ridge_prep 

10.  t2ugrid — I — bound 

I 

+-to_ugrid 

11.  abort_ice — I — ice_domain 

I 

H — ice_f ileunits 
I 

H — ice_mpi_internal 
I 

+- (shr_sys_abort) 

I 

H — end  run 


4.2.2  Directory  Structure 


The  present  PIPS  3.0  code  distribution  includes  makefiles,  input  files  and  scripts.  The  primary 
directory  is  pips3/,  and  a  run  directory  (rundir)  is  created  upon  the  initial  run  of  the  comp_ice 
script. 

+-pips3/  -  primary  directory 

+-README_V3.1  -  basic  information 
+-bld/  -  makefiles 

+-Macros.<OS>  -  macro  definitions  for  the  given  operating  system,  used  by 
Makefile.<OS> 

+-Makefile.<OS>  -  primary  makefile  for  a  given  operating  system 
(<std>  works  for  most  systems) 

+-makedep.c  -  perl  script  that  detennines  module  dependencies 
+-clean_ice  -  script  that  removes  files  from  the  compile  directory 
+-comp_ice  -  script  that  sets  up  the  run  directory  and  compiles  the  code 
+-doc/  -  documentation 

+-PIPS_SDD.doc-  software  design  description 
+-PIPS_UM.doc-  user’s  manual 
+-PIPS_VTR.doc-  validation  test  report 

+-cicedoc.pdf  -  CICE:  the  Los  Alamos  Sea  Ice  Model  Doc  and  User's  Manual 
+-PDF/  -  PDF  documents  of  several  publications  related  to  CICE 
+-ice.log.<OS>  -  sample  diagnostic  output  fdes 
+-source/  -  PIPS  3.0  source  code. 
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+-PIPS.F  -  main  program 
+-PIPS.F_debug-  debugging  version  of  PIPS. F 
+-ice_albedo.F  -  albedo  parameterization 

+-ice_atmo.F  -  stability-based  parameterization  for  calculation  of  turbulent  ice-atm  fluxes 

+-ice_calendar.F  -  keeps  track  of  what  time  it  is 

+-ice_constants.F  -  physical  and  numerical  constants  and  parameters 

+-ice_coupling.F  -  interface  with  the  flux  coupler 

+-ice_diagnostics.F  -  miscellaneous  diagnostic  and  debugging  routines 

+-ice_domain.F  -  MPI  subdomain  sizes  and  related  parallel  processing  info 

+-ice_dyn_evp.F  -  elastic-viscous-plastic  dynamics  component 

+-ice_exit.F  -  aborts  the  model,  printing  an  error  message 

+-ice_fileunits.F  -  unit  numbers  for  I/O 

+-ice_flux.F  -fluxes  needed/produced  by  the  model 

+-ice_flux_in.F  -  Reads  and  interpolates  forcing  data  for  stand-alone  ice  model  runs 
+-ice_grid.F  -  grid  and  land  masks 

+-ice_history.F  -  netCDF  output  routines  and  restart  read/write 
+-ice_init.F  -  namelist  and  initializations 
+-ice_itd.F  -  utilities  for  managing  ice  thickness  distribution 
+-ice_itd_linear.F  -  linear  remapping  for  transport  in  thickness  space 
+-ice_kinds_mod.F  -  basic  definitions  of  reals,  integers,  etc. 

+-ice_mechred.F  -  mechanical  redistribution  component  (ridging) 

+-ice_model_size.F  -  grid  size  and  number  of  thickness  categories  and  vertical  layers 

+-ice_model_size.F.gxl-  specific  ice_model_size.F  for  use  by  scripts  with  <  1°  >grid 

+-ice_model_size.Fgx3-  specific  ice_model_size.F  for  use  by  scripts  with  <  3°  >  grid 

+-ice_mpi_internal.F  -  utilities  for  internal  MPI  parallelization 

+-ice_ocean.F  -  mixed  layer  ocean  model 

+-ice_read_write.F  -  utilities  for  reading  and  writing  files 

+-ice_sca!ing.F  -  ice-area  scaling  of  variables  for  the  coupler 

+-ice_state.F  -  essential  arrays  to  describe  the  state  of  the  ice 

+-ice_therm_itd.F  -  thermodynamic  changes  related  to  ice  thickness  distribution  (post¬ 
coupling) 

+-ice_therm_vertical.F  -  vertical  growth  rates  and  fluxes  (pre-coupling  thennodynamics) 
+-ice_timers.F  -  timing  routines 

+-ice_transport_mpdata.F  -  horizontal  advection  via  MPDATA  or  upwind 
+-ice_transport_remap.F  -  horizontal  advection  via  incremental  remapping 
+-ice_work.F  -  globally  accessible  work  arrays 
+-rundir/  -  execution  or  "run"  directory  generated  when  the  code  is  compiled 
using  the  comp  ice  script 
+-cice  -  code  executable 

+-compile/  -  directory  containing  object  files,  etc. 

+-grid  -  horizontal  grid  file  from  pips/input_templates/ 

+-ice.log.[ID]  -  diagnostic  output  file 
+-ice_in  -  namelist  input  data  from  pips/input  templates 
+-hist/iceh_mavg.[timeID].nc  -  monthly  average  output  history  file 
+-kmt  -  land  mask  file  from  pips/input_templates/ 
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+-run_ice  -  batch  run  script  file  from  pips/input_templates/ 

5.0  PIPS  3.0  DETAILED  DESIGN 

The  following  sections  give  a  detailed  description  of  the  purpose,  variables,  logic,  and  constraints 

for  the  software  elements  in  the  model. 

5.1  Constraints  and  Limitations 

1 .  Fluxes  sent  to  the  coupler  could  have  incorrect  values  in  grid  cells  that  transform  from  an 
ice-free  state  to  having  ice  during  the  given  time  step,  or  vice  versa,  due  to  scaling  by  the 
ice  area.  The  flux  coupler  must  have  area  scaling  so  that  the  ice  and  land  models  are 
treated  reliably  in  the  coupler  (but  note  that  the  land  area  does  not  suddenly  become  zero 
in  a  grid  cell,  as  does  the  ice  area). 

2.  A  significant  fraction  (more  than  10%)  of  the  total  shortwave  radiation  is  absorbed  at  the 
surface.  It  should,  however,  be  penetrating  into  the  ice  interior  instead.  This  is  due  to  use 
of  the  aggregated,  effective  albedo  rather  than  the  bare  ice  albedo  when  snowpatch  < 

1.  Repairing  the  problem  will  require  more  albedo  arrays  to  be  added  to  the  code. 

3.  The  date-of-onset  diagnostic  variables,  melt_onset  and  f  rz_onset,  are  not 
included  in  the  restart  file.  These  could  therefore  be  incorrect  for  the  current  year  if  the 
run  is  restarted  after  Jan  1.  Also,  these  variables  were  created  with  the  Arctic  in  mind  and 
may  be  incorrect  for  the  Antarctic. 

4.  Timers  are  architecture  dependent. 

5.  Local  domains  are  not  padded  for  uneven  division  of  the  global  domain. 


5.2  Logic  and  Basic  Equations 

5.2.1  Coupling  with  Ocean  Model  Components 

PIPS  3.0  exchanges  information  with  the  NCOM  ocean  model  via  local  file  transfer.  The  fluxes 
and  state  variables  passed  between  the  PIPS  3.0  model  and  the  ocean  model  are:  ice 
concentration,  SST,  heat  flux  ice  to  ocean,  and  the  ocean/ice  stresses.  The  state  variables  passed 
between  the  ocean  model  and  PIPS  3.0  are:  ocean  SST  and  ocean  currents. 

By  convention,  directional  fluxes  are  positive  downward. 

The  ice  fraction  at  (aice ;  see  Section  6.0  for  a  description  of  all  typewritten  equivalents)  is  the 
total  fractional  ice  coverage  of  a  grid  cell.  That  is,  in  each  cell, 

ai  =  0  if  there  is  no  ice 

at  =1  if  there  is  no  open  water 

0  <  a,  <  1  if  there  are  both  ice  and  open  water 
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where  a,  is  the  sum  of  fractional  ice  areas  for  each  category  of  ice.  The  ice  fraction  is  employed 
by  the  flux  coupler  to  merge  fluxes  from  PIPS  3.0  with  fluxes  from  the  other  components.  For 
instance,  the  penetrating  shortwave  radiation  flux,  weighted  by  a,-,  is  combined  with  the  net 
shortwave  radiation  flux  through  ice-free  leads,  weighted  by  ( 1  -  a,),  to  get  the  net  shortwave 
flux  into  the  ocean  over  the  entire  grid  cell.  The  flux  coupler  requires  the  fluxes  to  be  divided  by 
the  total  ice  area  so  that  the  ice  and  land  models  are  managed  identically  (land  also  may  occupy 
less  than  100%  of  an  atmospheric  grid  cell).  These  fluxes  are  “per  unit  ice  area”  rather  than  “per 
unit  grid  cell  area.” 

5.2. 1.1  Atmosphere 

Wind  velocity,  air  density,  specific  humidity  and  potential  temperature  at  the  given  level  height 
zD  are  used  to  calculate  transfer  coefficients  in  formulas  for  the  surface  wind  stress  and  turbulent 
heat  fluxes  fa  ,Fs  and  Ft ,  as  described  below.  Wind  is  quite  possibly  the  main  forcing  mechanism 

for  the  ice  motion,  although  the  ice-ocean  stress,  Coriolis  force,  and  slope  of  the  ocean  surface 
are  also  important  [41].  The  sensible  and  latent  heat  fluxes,  Fs  and  Fi,  along  with  shortwave  and 
longwave  radiation,  F, ,  FjP  and  F  ^ ,  are  included  in  the  flux  balance  that  establishes  the  ice  or 

snow  surface  temperature.  As  described  in  Section  5. 2.2. 5,  these  fluxes  rely  nonlinearly  on  the 
ice  surface  temperature  Tsfc  The  balance  equation  is  iterated  until  convergence,  and  the  resulting 
fluxes  and  Tsfc  are  then  passed  to  the  flux  coupler.  The  flux  of  water  evaporated  to  the  atmosphere 
is  given  in  terms  of  the  latent  heat,  Feva p  =  FJ  (Lvap  +  Lice),  where  Lvap  and  Lice  are  latent  heats  of 
vaporization  and  fusion. 

The  snowfall  precipitation  rate  (specified  as  liquid  water  equivalent  and  converted  by  PIPS  3.0  to 
snow  depth)  also  contributes  to  the  heat  and  water  mass  budgets  of  the  ice  layer.  Although  melt 
ponds  usually  develop  on  the  ice  surface  in  the  Arctic  and  refreeze  later  in  the  autumn,  reducing 
the  total  amount  of  fresh  water  that  reaches  the  ocean  and  altering  the  heat  budget  of  the  ice, 
there  currently  is  no  melt  pond  parameterization.  Rain  and  all  melted  snow  end  up  in  the  ocean. 

Wind  stress  and  transfer  coefficients  for  the  turbulent  heat  fluxes  are  computed  in  subroutine 
atmo_boundary_layer  following  [23].  For  clarity,  the  equations  are  replicated  here  in  the  present 
notation.  The  wind  stress  and  turbulent  heat  flux  calculation  accounts  for  both  stable  and  unstable 
atmosphere-ice  boundary  layers.  Define  the  “stability”  as 

T_^.[  o'  ,  Q' 

,/J  0,(1  +  0.6068.)  1/0.606+8,; 

where  k  is  the  von  Kannan  constant,  g  is  gravitational  acceleration,  and  u*,  0*  and  Q*  are 
turbulent  scales  for  velocity,  temperature  and  humidity,  respectively: 
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=  ce(®.~T4c)  0) 

Q‘  =  £;,  ( Q,  -  £?* )  • 

where  the  wind  speed  has  a  minimum  value  of  1  m/s.  Ice  motion  is  ignored  in  u  ,  and  r*  and 
Qsfc  represent  the  surface  temperature  and  specific  humidity,  respectively.  The  latter  is  calculated 
by  assuming  a  saturated  surface  temperature  Tsfc,  as  described  in  Section  5.2.2.5.I. 

The  exchange  coefficients  cu,  cq  and  cq  are  initialized  as 

K 

Hzref/zice) 

and  updated  during  a  short  iteration,  as  they  depend  on  the  turbulent  scales.  Here,  zref  is  a 
reference  height  of  10  m  and  zice  represents  the  roughness  length  scale  for  the  given  sea  ice 

0  25 

category.  T is  constrained  to  have  a  magnitude  less  than  10.  Further,  defining  %  =  ( 1  -  16T)  ' 
and  y  >1 ,  the  “integrated  flux  profiles”  for  momentum  and  stability  in  the  unstable  ( T  <  0)  case 
are  given  by 

Vm  =  21n[0.5(l  + j)]  +  ln[0.5(l  +  r)]-2tan-1J  +  | 

¥s=  21n[0.5(l  +  r)]. 

Aside  from  the  parameterization  used  in  [23],  profiles  are  used  for  the  stable  case  from  [22], 

¥,„  =¥,  =-[0.7T  +  0.75  (T- 14.3)  exp  (-0.35T)  + 10.7], 

The  coefficients  are  then  updated  as 

cu 

1  +  cu(X-y/m)lK 

_ £/? _ 

1  +  cg(X-y/s)lK 


where  X  =  In (z0  /zre f ).  The  first  iteration  is  completed  with  new  turbulent  scales  from  equation 
(1).  After  five  iterations  the  latent  and  sensible  heat  flux  coefficients  are  computed,  along  with 
the  wind  stress: 
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C,  =  Pa(LVap+Lice)U*Cq 

Cs  =  pacpuce+ 1, 


where  pa  is  the  density  of  air  and  cp  is  its  specific  heat.  Using  [22]  again,  a  constant  has  been 
added  to  the  sensible  heat  flux  coefficient  to  allow  some  heat  to  pass  between  the  atmosphere  and 
the  ice  surface  in  calm,  stable  conditions.  The  atmospheric  reference  temperature  Tp1  is 

computed  from  Ta  and  Ts/C  using  the  coefficients  cu,  cq  and  cq.  Although  PIPS  3.0  does  not  use 
this  quantity,  it  is  quite  convenient  for  the  ice  model  to  execute  this  calculation.  The  atmospheric 
reference  temperature  is  returned  to  the  flux  coupler  as  a  climate  diagnostic.  The  same  holds  true 
for  the  reference  humidity, 

Additional  details  about  the  latent  and  sensible  heat  fluxes  and  other  quantities  referred  to  here 
can  be  found  in  Section  5. 2.2. 5.1. 


5. 2. 1.2  Ocean 

New  sea  ice  grows  when  the  ocean  temperature  drops  below  its  freezing  temperature,  7>  =  -pS, 
where  S  is  the  seawater  salinity  and  p  =0.054  is  the  ratio  of  the  freezing  temperature  of  brine  to 
its  salinity.  The  ocean  models,  either  NCOM  or  HYCOM,  perform  this  calculation;  if  the 
freezing/melting  potential  Ffrzm[t  is  positive;  its  value  is  a  certain  quantity  of  frazil  ice  that  has 
formed  in  one  or  more  layers  of  the  ocean  and  floated  to  the  surface.  (The  ocean  models  assume 
that  the  quantity  of  new  ice  implied  by  the  freezing  potential  actually  forms.)  Generally,  this  ice 
is  added  to  the  thinnest  ice  category.  The  new  ice  is  created  in  the  open  water  area  of  the  grid  cell 
to  a  specified  minimum  thickness.  Therefore,  if  the  open  water  area  is  nearly  zero  or  if  there  is 
more  new  ice  than  will  fit  into  the  thinnest  ice  category,  the  new  ice  is  then  spread  unifonnly 
over  the  entire  grid  cell. 

If  Ffrzmit  is  negative,  it  is  used  to  heat  existing  ice  from  below.  Specifically,  the  sea  surface 
temperature  and  salinity  are  used  to  compute  an  oceanic  heat  flux  Fw  ( |  Fw  |  <  |  Ffrzmlt  | )  which  is 

applied  at  the  bottom  of  the  ice.  The  fraction  of  the  melting  potential  actually  used  to  melt  ice  is 
returned  to  the  coupler  in  F]met.  The  ocean  models  adjust  their  own  heat  budgets  with  this 
quantity,  assuming  that  the  rest  of  the  flux  stayed  in  the  ocean. 

In  addition  to  runoff  from  rain  and  melted  snow,  the  fresh  water  flux  Fwater  includes  ice 
meltwater  from  the  top  surface  of  the  ice  and  water  melted  or  frozen  (a  negative  flux)  at  the 
bottom  surface.  This  flux  is  computed  as  the  net  change  of  fresh  water  in  the  ice  and  snow 
volume  over  the  coupling  time  step,  prohibiting  frazil  ice  fonnation  and  newly  accumulated 
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snow. 

There  is  a  flux  of  salt  into  the  sea  under  melting  conditions,  and  a  (negative)  flux  when  ocean 
water  is  freezing.  However,  melting  sea  ice  ultimately  freshens  the  top  layer  of  the  ocean,  since  it 
is  much  more  saline  than  the  ice.  The  PIPS  3.0  model  passes  the  net  flux  of  salt,  Fsait,  to  the  flux 
coupler,  based  on  the  net  exchange  in  salt  for  ice  in  every  category.  In  the  present  configuration, 
ice_ref_salinity  is  used  in  computing  the  salt  flux,  although  the  ice  salinity  utilized  in  the 
thennodynamic  calculation  has  differing  values  in  the  ice  layers. 

A  fraction  of  the  incoming  shortwave  F  n  penetrates  the  snow  and  ice  layers  and  passes  into  the 
ocean,  as  described  in  Section  5.2.2.5.I. 


Many  ice  models  compute  the  sea  surface  slope  VHo  from  geostrophic  ocean  currents  supplied 
by  an  ocean  model  or  other  data  source.  In  this  case,  the  sea  surface  height  Ho  is  a  prognostic 

variable  in  the  NCOM  and  HYCOM  models.  The  flux  coupler  provides  the  surface  slope 
directly,  instead  of  inferring  it  from  the  currents  (The  option  of  computing  surface  slope  from  the 
currents  is  provided  in  subroutine  evp _prep.).  The  PIPS  3.0  model  uses  the  surface  layer  currents 
Uw  to  derive  the  stress  between  the  ocean  and  the  ice,  and  subsequently  the  ice  velocity  u  .  This 
stress,  in  relation  to  the  ice, 

{Uw  —  u\ cos 6  +  k x {Uw  - u ) sin 6 

is  then  transferred  to  the  flux  coupler  (relative  to  the  ocean)  for  use  by  the  ocean  model.  Here,  6 
is  the  turning  angle  between  geostrophic  and  surface  currents,  cw  is  the  ocean  drag  coefficient, 

pw  is  the  density  of  seawater  (dragw=  cwpw),  and  k  represents  the  vertical  unit  vector.  The 

turning  angle  is  necessary  if  the  top  ocean  model  layers  cannot  resolve  the  F.kman  spiral  in  the 
boundary  layer.  If  the  top  layer  is  sufficiently  thin  compared  to  the  typical  depth  of  the  Ekman 
spiral,  then  6  =0  is  a  good  approximation.  Here  it  is  assumed  that  the  top  layer  is  sufficiently 
thin. 


CwP» 


U,-u 


5.2.2  Model  Components 

The  Arctic  and  Antarctic  sea  ice  packs  are  comprised  of  combinations  of  open  water,  thin  first- 
year  ice,  thicker  multi-year  ice,  and  thick  pressure  ridges.  The  thennodynamic  and  dynamic 
characteristics  of  the  ice  pack  depend  on  the  amount  of  ice  lying  in  each  thickness  range.  Thus 
the  basic  challenge  in  sea  ice  modeling  is  to  describe  the  evolution  of  the  ice  thickness 
distribution  (ITD)  in  both  time  and  space. 


The  fundamental  equation  solved  by  PIPS  3.0  is  from  [43]: 


37 


PIPS  3.0  SDD 


dg  d 

-£  =  -V.(gu)--(fg)  +  y,, 

ot  oh 


(2) 


where  u  is  the  horizontal  ice  velocity,  V  =  ( — ,  — )  ,  /  is  the  rate  of  thermodynamic  ice  growth, 

dx  dy 

\fj  is  a  ridging  redistribution  function,  and  g  represents  the  ice  thickness  distribution  function. 
We  define  g(x,h,t)dh  as  the  fractional  area  covered  by  ice  in  the  thickness  range  ( h,h  +  dh )  at  a 
given  time  and  location. 

Equation  (2)  is  solved  by  sectioning  the  ice  pack  at  each  grid  point  into  discrete  thickness 
categories.  The  number  of  categories  is  set  by  the  user,  with  a  default  value  Nc  =5.  (Five 
categories,  plus  open  water,  are  sufficient  to  simulate  the  yearly  cycles  of  ice  thickness,  ice 
strength,  and  surface  fluxes  [6],  [25].  Each  category  n  has  lower  thickness  bound  Hn_ i  and 

upper  bound  Hn.  The  lower  bound  of  the  thinnest  ice  category,  H0,  is  set  to  zero.  The  other 
boundaries  are  selected  with  greater  resolution  for  small  h,  since  the  properties  of  the  ice  pack 
are  particularly  sensitive  to  the  amount  of  thin  ice  [28].  The  continuous  function  g(h)  is 
replaced  by  the  discrete  variable  ain,  defined  as  the  fractional  area  covered  by  ice  in  the  thickness 
range  ( Hn  [,Hn).  The  fractional  area  of  open  water  is  denoted  by aj0,  giving  ain  =  1  by 
definition. 

Category  boundaries  are  calculated  in  initjtd  using  one  of  two  formulas.  The  old  fonnula,  from 
[26],  assigns  lower  boundaries  (in  meters)  of  (0.0,  0.64,  1.39,  2.47,  and  4.57)  for  categories  1  to 
5  when  N_C  =  5.  A  new  formula  has  been  created  for  boundaries  that  are  round  numbers.  This 
formula  gives  boundaries  (0.0,  0.60,  1.40,  2.40,  and  3.60)  for  N_C  =  5.  The  old  formula 
(kcatbound  =  0  in  the  namelist)  is  the  default.  A  user  may  substitute  his  or  her  own  preferred 
boundaries  in  init  itd. 


Besides  the  fractional  ice  area,  ain,  the  following  state  variables  for  each  category  n  are  defined 
as: 

•  vin,  the  ice  volume,  equal  to  the  product  of  ain  and  the  ice  thickness  hin. 

•  vsn,  the  snow  volume,  equal  to  the  product  of  ain  and  the  snow  thickness  hsn. 

•  Oink,  the  internal  ice  energy  in  layer  k,  equal  to  the  product  of  the  ice  layer  volume, 

and  the  ice  layer  enthalpy,  qink.  Here  A)  is  the  total  number  of  ice  layers,  with  a  default 
value  Nj  =4,  and  qink  is  the  negative  of  the  energy  needed  to  melt  a  unit  volume  of  ice  and 
raise  its  temperature  to  0  C.  This  is  discussed  in  Section  5. 2.2. 5.  (NOTE:  In  the  current 
code,  et  <  0  and  qt  <  0  with  ej  =  viqj .) 

•  esn,  the  snow  energy,  equal  to  the  product  of  vsn  and  the  snow  enthalpy,  qsn.  At  this 


38 


PIPS  3.0  SDD 


writing  there  is  only  one  snow  layer,  but  future  versions  of  PIPS  3.0  will  allow  for 
multiple  snow  layers.  (Similarly,  es  <  0  in  the  code.) 

•  Tsfn,  the  surface  temperature. 

3 

Because  the  fractional  area  is  unitless,  the  volume  variables  have  units  of  meters  (i.e.,  m  of  ice 
2  2 
or  snow  per  m  of  grid  cell  area),  and  the  energy  variables  have  units  of  J/m  . 

The  three  tenns  on  the  right-hand  side  of  equation  (2)  illustrate  three  kinds  of  sea  ice  transport: 

(1)  horizontal  transport  in  (x,y)  space; 

(2)  transport  in  thickness  space  h  due  to  thermodynamic  growth  and  melting;  and 

(3)  transport  in  thickness  space  h  due  to  ridging  and  other  mechanical  processes. 

The  equation  is  solved  by  operator  splitting  in  three  stages,  with  two  out  of  three  terms  on  the 
right  set  to  zero  in  each  stage.  Horizontal  transport  is  computed  using  the  incremental  remapping 
scheme  of  [9]  as  adapted  for  sea  ice  by  [26].  This  scheme  is  discussed  in  Section  5.2.2. 1.  Ice  is 
carried  in  thickness  space  through  the  remapping  scheme  of  [25],  as  described  in  Section  5.2. 2.2. 
The  mechanical  redistribution  scheme,  based  on  [43],  [35],  [16],  [13],  and  [27]  is  outlined  in 
Section  5. 2.2. 3.  To  solve  the  horizontal  transport  and  ridging  equations,  we  need  the  ice  velocity 
u.  To  calculate  transport  in  thickness  space,  ice  growth  rate  /  must  be  known  in  each  thickness 
category.  The  elastic-viscous-plastic  (EVP)  ice  dynamics  scheme  of  [18]  is  used,  as  modified  by 
[17],  [19]  and  [20]  to  find  the  velocity,  described  in  Section  5. 2.2. 4.  Finally,  the  thermodynamic 
model  of  [7],  discussed  in  Section  5.2. 2. 5,  is  used  to  compute  f. 


5.2. 2.1  Horizontal  Transport 

Solving  the  continuity  or  transport  equation, 

%-+V-(a„a)  =  0,  (3) 

at 

for  the  fractional  ice  area  in  each  thickness  category  n.  Equation  (3)  describes  the  conservation  of 
ice  area  in  horizontal  transport.  It  is  taken  from  Equation  (2)  by  discretizing  g  and  ignoring  the 
second  and  third  terms  on  the  right-hand  side,  which  are  handled  separately  (see  Sections  5. 2. 2. 2 
and  5. 2.2. 3). 

There  are  comparable  conservation  equations  for  ice  volume,  snow  volume,  ice  energy,  snow 
energy,  and  area-weighted  surface  temperature: 
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— TT"  "I-  V '  (V/„u)  =  0, 
ot 

(4) 

s"  +V-(vOTu)  -  0, 

Ot 

(5) 

r)p 

-#  +  V-(eMu)  =  0, 

Ot 

(6) 

8c 

"  +V-(es„u)  -  0, 
ot 

(V) 

8(^)+v  =  „ 

(8) 

For  simplicity,  ice  and  snow  are  assumed  to  have  constant  densities,  so  that  volume  conservation 
is  equal  to  mass  conservation.  Also  transported  is  the  fractional  area  of  open  water,  using 
equation  (3)  with  n  =0.  Including  Nc=  5  and  Nt  =  4  there  are  46  transport  equations  to  be 
solved. 

Three  transport  schemes  are  available,  upwind,  MPDATA  [39]  and  the  incremental  remapping 
scheme  of  [9]  as  modified  for  sea  ice  by  [26].  Because  several  fields  are  transported,  the  transport 
module  is  quite  computationally  expensive  (almost  half  the  total  computer  time)  in  runs  using 
MPDATA.  Although  a  cheaper  first-order  upwind  scheme  is  available  as  an  MPDATA  option 
(see  Section  4. 1.3.1),  it  is  advised  that  the  incremental  remapping  method  be  used  instead.  This 
scheme  has  many  desirable  features: 

•  It  conserves  the  quantity  being  transported  (area,  volume,  or  energy). 

•  It  is  non-oscillatory,  meaning  it  does  not  create  spurious  ripples  in  the  transported  fields. 

•  It  preserves  tracer  monotonicity,  meaning  it  does  not  create  new  extrema  in  the  thickness 
and  enthalpy  fields.  The  values  at  time  m  +1  are  bounded  by  the  values  at  time  m. 

•  It  is  second-order  accurate  in  space  and  therefore  is  a  great  deal  less  diffusive  than  first- 
order  schemes.  The  accuracy  can  be  decreased  locally  to  first  order  to  preserve 
monotonicity. 

•  It  is  efficient  for  large  amounts  of  categories  or  tracers.  Much  of  the  work  is  geometrical 
and  is  executed  only  once  per  grid  cell  instead  of  being  repeated  for  each  quantity  being 
transported. 

The  time  step  is  limited  by  the  requirement  that  trajectories  projected  backward  from  grid  cell 
comers  are  confined  to  the  four  surrounding  cells,  thus  defining  incremental  remapping  as 
opposed  to  general  remapping.  This  requirement  leads  to  a  CFL-like  condition, 

max  |  u  |  At 
Ax 
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For  greatly  divergent  velocity  fields  the  maximum  time  step  must  be  decreased  by  a  factor  of  two 
to  ensure  that  trajectories  do  not  cross.  Then  again,  ice  velocity  fields  in  climate  models  typically 
have  small  divergences  per  time  step  relative  to  the  grid  size. 

The  remapping  algorithm  may  be  summarized  as  follows: 

1.  Given  mean  values  of  the  ice  area  and  tracer  fields  in  each  grid  cell,  generate  linear 
approximations  of  these  fields.  Limit  the  field  gradients  to  conserve  monotonicity. 

2.  Given  ice  velocities  at  grid  cell  corners,  identify  departure  regions  for  the  fluxes  across 
each  cell  edge.  Divide  these  departure  regions  into  triangles  and  compute  the  coordinates 
of  the  triangle  vertices. 

3.  Integrate  these  fields  across  the  departure  triangles  to  obtain  the  area,  volume,  and  energy 
fluxes  across  each  cell  edge. 

4.  Transfer  the  fluxes  over  cell  edges  and  update  the  state  variables. 

Because  all  scalar  fields  are  transported  by  the  same  velocity  field,  step  (2)  is  performed  only 
once  per  time  step.  The  other  three  steps  are  repeated  for  each  field  in  each  thickness  category. 
These  steps  are  described  below. 


5.2.2. 1.1  Reconstructing  Area  and  Tracer  Fields 

Firstly,  using  the  known  values  of  the  state  variables,  the  ice  area  and  tracer  fields  are 
reconstructed  in  every  grid  cell  as  linear  functions  of  x  and  y.  For  each  field  the  value  at  the  cell 
center  is  computed  (i.e.,  at  the  origin  of  a  2D  Cartesian  coordinate  system  defined  for  that  grid 
cell),  along  with  gradients  in  the  x  and  y  directions.  The  gradients  are  limited  to  conserve 
monotonicity.  When  integrated  over  a  grid  cell,  the  reconstructed  fields  must  have  mean  values 

equal  to  the  known  state  variables,  given  by  a  for  fractional  area,  h  for  thickness,  and  q  for 
enthalpy.  The  mean  values  are  typically  not  equal  to  the  values  at  the  cell  center.  For  instance, 
the  mean  ice  area  must  equal  the  value  at  the  centroid,  which  may  not  be  at  the  cell  center. 

First  consider  the  fractional  ice  area,  the  analog  to  fluid  density  p  in  [9],  For  each  thickness 
category  a  field  a(r)  is  constructed  whose  mean  is  a  ,  where  r  =(x,  y )  is  the  position  vector 
relative  to  the  cell  center.  That  is,  we  require 

[  adA  =  aA,  (9) 

J  A 

where  A  =  [  dA  is  the  grid  cell  area.  Equation  (9)  is  satisfied  if  a(r)  takes  the  form 


a(r)  =  a  +aa  <  Va  >  -(r-r),  (10) 

where  <Va>  is  the  centered  estimate  of  the  area  gradient  in  the  cell,  aa  is  a  limiting  coefficient 

that  implements  monotonicity,  and  r  is  the  cell  centroid: 
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r  =  —  [  rdA. 

A3  a 

It  follows  from  (10)  that  the  ice  area  at  the  cell  center  (r  =0)  is 

ac  =a-axx-ayy. 


where  ax  =  aa(da/Sx)  and  a  =  ajca/dy)  are  the  limited  gradients  in  the  x  and  y  directions, 
respectively,  and  the  components  of  r,  x  =  jxdA/A  and  v  =  j  vdA/A  are  tested  using  the 
triangle  integration  formulas  described  in  Section  5. 2.2. 1.3.  These  means,  combined  with  higher 
order  means  such  as  x2 ,  xy,  and  y2 ,  are  computed  once  and  stored. 


Next  consider  the  snow  and  ice  thickness  and  enthalpy  fields.  Thickness  is  analogous  to  the 
tracer  concentration  T  in  [9]  but  there  is  no  analog  in  their  study  to  the  enthalpy.  The 
reconstructed  ice  or  snow  thickness  h(r)  and  enthalpy  q( r)  must  satisfy 


ahdA  =  ah  A, 

J  A 

ahqdA  =  ahqA. 

J  A 


(ID 

(12) 


Equations  (11)  and  (12)  are  satisfied  when  h( r)  and  q( r)  are  denoted  by 


h( r)  =  h  +  ah  <Vh>  ■( r  - r ), 


(13) 


q(r)  =  q  +  a<Vq>{r-r), 


(14) 


where  ah  and  a  are  limiting  coefficients,  r  is  the  center  of  ice  area, 


—  f  ardA, 


a  A  Ja 


and  r  is  the  center  of  ice  or  snow  volume, 


r  =  — tL-  f  ahrdA. 
ahA  3 A 


Evaluating  the  integrals,  we  find  that  the  components  of  r  are  given  by 
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acx  +  axx 2  +a  xy 

K  = - — - - - , 

a 

~  acy  +  axxy  +  a  y1 

y  = - = — - — 

a 


and  the  components  of  r  are 


where 


clx  +  c2x 2  +  c3xy +  c4x3  +  c5x2y  +  c6xy2 


>  _cxy  +  c2 xy  +  c3  y 2  +  c4 x2y  +  c5 xy2  +c6y3 
y  cih 


c, 

c2 

C3 

C4 


acK, 

achx+ajic, 

achy+ayhc. 


aA+a>K 

aA- 


From  equations  (13)  and  (14),  the  thickness  and  enthalpy  at  the  cell  center  are  denoted  by 

h  =  h  -hxx-hry, 

qc  =  q-qxx-qyy. 


where  hx,  hy,  qx  and  qv  are  the  limited  gradients  of  thickness  and  enthalpy.  The  surface 
temperature  is  dealt  with  the  same  way  as  ice  or  snow  thickness,  but  it  has  no  associated 
enthalpy. 

Monotonicity  is  preserved  by  van  Leer  limiting.  If  <j>(i,j)  is  the  mean  value  of  some  field  in 
grid  cell  (/,  /),  centered  gradients  of  (f)  in  the  x  and  y  directions  are  computed.  They  are  also 
checked  to  see  that  the  gradients  provide  values  of  (f>  within  cell  (/,  /)  that  lie  outside  the  range  of 
(f>  in  the  grid  cell  and  its  eight  neighbors.  Let  (f)mix  and  <f)mm  be  the  maximum  and  minimum 
values  of  (f>  over  the  cell  and  its  neighbors,  and  let  (f)m.ix  and  (f>mm  be  the  maximum  and  minimum 
values  of  the  reconstructed  6  within  the  cell.  Since  the  reconstruction  is  linear,  (b  and  (f> . are 
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located  at  cell  comers.  If  (j)mm  >  ^max  or  <f>min  <  <j)mm ,  the  unlimited  gradient  is  multiplied  by 
a  =  min(«max ,  amm  )  ,  where 

«max  =(4ax-^y(^max-A 

a  =  (6  -(j))K(j)  .  -<b). 

min  Vrmm  V  /'  \  Y  niin  V  /* 

Otherwise  the  gradient  does  not  need  to  be  limited. 


5.2. 2. 1.2  Locating  Departure  Triangles 

The  technique  for  locating  departure  triangles  is  discussed  in  detail  by  [9].  The  basic  idea  is  seen 
in  Figure  1,  which  illustrates  a  shaded  quadrilateral  departure  region  whose  contents  are  carried 
to  the  target  or  home  grid  cell,  labeled  H.  The  neighboring  grid  cells  are  tagged  by  compass 
directions:  NW,  N,  NE,  W,  and  E.  The  four  vectors  point  along  the  velocity  field  at  the  cell 
corners,  and  the  departure  region  is  formed  by  joining  the  starting  points  of  the  vectors.  Rather 
than  integrating  across  the  entire  departure  region,  we  compute  fluxes  across  cell  edges. 
Departure  regions  are  selected  for  the  north  and  east  edges  of  each  cell,  which  are  also  the  south 
and  west  edges  of  neighboring  cells.  Consider  the  north  edge  of  the  home  cell,  over  which  there 
are  fluxes  from  the  neighboring  NW  and  N  cells.  The  contributing  region  from  the  NW  cell  is  a 
triangle  with  vertices  abc,  and  that  region  from  the  N  cell  is  a  quadrilateral  that  can  be  divided 
into  two  triangles  with  vertices  acd  and  ade.  Focusing  on  triangle  abc,  the  coordinates  of  vertices 
b  and  c  with  respect  to  to  the  cell  corner  (vertex  a)  are  detennined  first,  using  Euclidean 
geometry,  to  find  vertex  c.  Then  these  three  vertices  are  translated  to  a  coordinate  system  located 
in  the  NW  cell  center.  This  translation  is  necessary  for  integrating  fields  (Section  5.2. 2. 1.3)  in  the 
coordinate  system  where  they  have  been  reconstructed  (Section  5. 2.2. 1.1).  Repeating  this  process 
for  the  north  and  east  edges  of  each  grid  cell,  we  compute  the  vertices  of  all  the  departure 
triangles  associated  with  each  cell  edge. 
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Figure  1:  In  incremental  remapping,  conserved  quantities  are  remapped  from  the  shaded 
departure  region,  which  is  a  quadrilateral  formed  by  joining  the  backward  trajectories 
from  the  four  cell  corners,  to  grid  cell  H.  The  region  fluxed  over  the  north  edge  of  cell  H 
consists  of  a  triangle  ( abc )  in  the  NW  cell  and  a  quadrilateral  (two  triangles,  acd  and  ade)  in 
the  N  cell. 

Figure  2,  reproduced  from  [9],  shows  all  possible  triangles  that  can  contribute  fluxes  across  the 
north  edge  of  a  grid  cell.  There  are  20  triangles,  which  may  be  divided  into  five  groups  of  four 
mutually  exclusive  triangles,  as  shown  in  Table  2.  In  this  table,  (xi,yi)  and  {xi,yi)  are  the 
Cartesian  coordinates  of  the  departure  points  relative  to  the  NW  and  NE  cell  corners, 
respectively.  The  departure  points  are  connected  by  a  straight  line  that  intersects  the  west  edge  at 
(0, ya)  relative  to  the  NW  comer  and  intersects  the  east  edge  at  (0,^) relative  to  the  NE  corner. 

This  line  intersects  the  N  edge  at  (xa,0)  relative  to  the  NW  corner  and  (xA,0)  relative  to  the  NE 
corner. 
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Figure  2:  The  20  possible  triangles  that  can  contribute  fluxes  across  the  north  edge  of  a  grid 
cell. 


From  Cartesian  geometry  it  can  be  shown  that 


yb 

x/, 


yAx+(x2y1-x1y2) 
Ax  +  x2  -  xl 
y2Ax  +  (x2yx  -  xt  v2 ) 
Ax +  x2-  Xj 

ya-yb  ’ 
yb*x 

ya~yb  ’ 


where  Ax  represents  the  length  of  the  N  edge.  The  east  cell  triangles  and  selecting  conditions  are 
alike  except  for  a  rotation  through  90  degrees. 
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Triangle  Group 

Triangle  Label 

Selecting  Logical  Condition 

1 

NW 

ya  >  0  and  y,  >  0  and  x,  <  0 

NW1 

ya<  0  and  y,  >  0  and  x,  <  0 

W 

ya  <  0  and  y,  <  0  and  x,  <  0 

W2 

ya  >  0  and  y,  <  0  and  x,  <  0 

2 

NE 

yh  >  0  and  y2  >  0  and  x2  >  0 

NE1 

yh  <  0  and  y2  >  0  and  x2  >  0 

E 

yh  <  0  and  y2  <  0  and  x2  >  0 

E2 

yh  >  0  and  y2  <  0  and  x2  >  0 

3 

W1 

ya<  0  and  y,  >  0  and  x,  <  0 

NW2 

ya  >  0  and  y,  <  0  and  x,  <  0 

El 

yh  <  0  and  y2  >  0  and  x2  >  0 

NE2 

yh  >  0  and  y2  <  0  and  x2  >  0 

4 

Hla 

yayb  - 0  and  >/ + >/,  < 0 

Nla 

yayb  -  0  and  ya  +  yb>  0 

Hlb 

yayb  <  ° and  y,  <  ° 

Nib 

yayb  <  0  and  y,  >  0 

5 

H2a 

yayb  -  ° and  >/ + >/,  < 0 

N2a 

yayb  -  0  and  ya  +  yb>  o 

H2b 

yayb  <  ° and  y2  <  ° 

N2b 

yayb  <  0  and  y2  >  0 

Table  4:  Evaluation  of  contributions  from  the  20  triangles  across  the  north  cell  edge.  The 
coordinates  x, ,  x2 ,  y, ,  y2 ,  ya ,  and  yb  are  defined  in  the  text.  We  define  y,  =  y1  if  Xi  >0, 
else  y,  =  ya .  Similarly,  y2  =  y2  if  x2  <  0  else  y2  =  yh . 

This  scheme  was  created  for  rectangular  grids.  Grid  cells  in  PIPS  3.0  essentially  lie  on  the 
surface  of  a  sphere  and  must  be  projected  onto  a  plane.  Many  projections  are  possible.  The 
projection  used  in  PIPS  3.0,  illustrated  in  Figure  3,  approximates  spherical  grid  cells  as 
quadrilaterals  in  the  plane  tangent  to  the  sphere  at  a  point  inside  the  cell.  The  quadrilateral 
vertices  are  {N/2,  E/2),  {-N/2,  W/2),  (-5/2,  -W/ 2),  and  (5/2,  -E/2),  where  N,  S,  E,  and  W  are  the 
lengths  of  the  cell  edges  on  the  spherical  grid.  The  quadrilateral  area,  (N  +  S){E  +  W  )/4,  is  an 
apt  approximation  to  the  true  spherical  area.  But  cell  edges  in  this  projection  are  not  orthogonal 
(i.e.,  they  do  not  converge  at  right  angles)  as  they  do  on  the  spherical  grid.  Therefore,  when 
vectors  are  translated  from  cell  comers  to  cell  centers,  the  departure  points  in  the  cell-center 
coordinate  system  must  be  inside  the  grid  cell  contributing  the  flux.  Otherwise,  monotonicity 
may  be  violated,  because  van  Leer  limiting  does  not  apply  outside  the  grid  cell. 
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f  j 


Figure  3:  A  grid  cell  on  the  surface  of  a  sphere  with  unequal  sides  of  length  N,  S,  E,  and  W 
is  estimated  as  a  quadrilateral  lying  in  the  tangent  plane  at  the  cell  center.  The 
quadrilateral  vertices  are  (N/2,  E/2),  (-N/2,  W/2),  (-S/2,  -W/2),  and  (S/2,  -E/2).  The  basis 
vectors  ( /  ,/ ),  located  at  the  northeast  cell  corner,  have  been  projected  into  the  cell-center 
coordinate  system  and  vary  from  the  cell-center  basis  vectors  ( i,j  ).  The  angles  6 N and  0E 
relating  the  two  bases  are  defined  in  the  text. 


Figure  3  demonstrates  the  difficulty.  At  the  cell  center  orthogonal  basis  vectors  i  and  j  that 
point  to  the  midpoints  of  the  cell  edges  are  defined.  Similarly,  at  each  cell  corner  is  defined  a 
coordinate  system  whose  basis  vectors,  i  and  j  point  along  cell  edges.  The  vectors  i  and  j 
are  orthogonal  in  the  cell-corner  reference  frame,  but  not  when  projected  into  the  reference  frame 
of  the  bordering  cell  center.  Because  of  this,  a  simple  transformation  is  used  to  preserve 

monotonicity  when  vectors  are  translated  from  corners  to  centers.  Consider  a  vector  (x  i  +  y  j  ) 
in  the  cell-comer  basis.  An  assumption  is  made  that  this  vector  has  the  same  coordinates  when  i 
and  j  are  non-orthogonal  projections  of  the  cell-corner  basis  vectors  into  the  cell-center  tangent 
plane,  as  in  Figure  3.  Thus  a  transfonnation  follows  from  the  ( i  ,j  )  basis  to  the  ( i ,  j  )  basis.  In 
the  cell-center  coordinate  system,  i  is  obtained  by  a  rotation  of  i  through  an  angle  0N ,  where 


0N  =  arctan 


(E-W \ 
2  N 


(15) 


Similarly,  j  is  obtained  by  a  rotation  of  j  through  0E ,  where 
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(  S-N \ 

6r=  arctan  -  .  (16) 

l  2  E  ) 

Vectors  are  transformed  between  basis  sets  using 

x)  f  cos  0N  -sin0£V  x 

y)  v  sin^v  cos^JU' 

which  may  be  verified  by  inspection,  alternately  setting  x  =  0  and  y  =  0 .  Similar 
transfonnations  are  used  at  the  three  other  cell  comers.  These  transfonnations  guarantee  that  the 
grid  cell  in  which  a  given  departure  point  is  positioned  does  not  change  under  a  modification  in 
coordinate  systems. 

Most  grids  cells  are  fairly  rectangular,  unlike  the  distorted  cell  shown  in  Figure  3.  On  the  1° 
displaced-pole  grid  frequently  used  for  PIPS  3.0  runs,  the  maximum  angle  in  equations  (15)  and 
(16)  is  about  1°.  Vector  transfonnations  could  therefore  be  left  out  on  most  grids  with  little  loss 
of  accuracy.  They  have  been  retained,  however,  because  they  assure  precise  monotonicity  at  little 
added  cost. 

Another  change  should  be  noted  in  the  scheme  of  [9]  for  locating  triangles.  In  that  paper, 
departure  points  are  defined  by  projecting  cell  comer  velocities  directly  backward.  That  is, 

xD  =  -u  A?,  (18) 

where  \D  is  the  location  of  the  departure  point  relative  to  the  cell  corner.  The  primes  represent 

vectors  defined  in  the  cell-comer  basis.  This  approximation  is  accurate  only  for  first-order.  The 
accuracy  may  be  raised  to  second-order  by  rectifying  the  velocity  with  a  midpoint  approximation 
before  finding  the  departure  point. 

First,  the  midpoint  of  the  backward  trajectory  is  estimated  using  xM  =  xD/2 ,  then  an  equation  like 
equation  (17)  is  used  to  transform  xM  to  the  appropriate  cell-center  coordinate  system.  Next,  a 
bilinear  interpolation  is  employed  to  estimate  the  velocity  at  xM  .  In  a  square  2x2  grid  cell  with 
cornersXj  =  (-1,-1)  ,x2  =  (1,-1),  x3  =  (l,l),and  x4  =  (-1,1),  the  values  of  any  scalar  field  <f) 
can  be  paired  up  at  the  cell  corners  with  the  following  bilinear  approximation: 

v  y)  =  - 1  )(y  - 1 )  -  fa  (* + 1  )(y  - 1 ) + k  (* + 1  )(y + 1 )  -  A  (*  - 1  )(y + 1  )L  (i  9) 


49 


PIPS  3.0  SDD 


where  (f\ ,  fa,  fa,  and  fa  are  the  corner  values.  To  use  equation  (19),  xMmust  be  transformed 
from  cell-center  coordinates  (x,y)  into  the  simpler  (x,  y)  coordinate  system.  Substituting  x  and 
y  for  (f>  in  equation  (19),  we  get 


x(x,  y)  =  ^-[x,  (x  - 1  )(y  - 1 )  -  x2  (x  + 1  )(y  - 1 )  +  x3  (x  + 1  )(y  + 1 )  -  x4  (x  - 1  )(y  + 1 )] , 
y(x,y)  =  I[Vi(x-l)(v-l)-y2(x  +  l)(y-l)  +  y3(T  +  l)(v  +  l)-y4(x-l)(y  +  l)]. 


Recalling  that  the  corner  coordinates  are  x,  =  {-S/2, -W  12.) ,  x2  =  (S/2, -E/2)  x3  =  (N/2,  E/2) , 
and  x4  =  (-N/2,  W/2) ,  expressions  are  derived  for  x  and  y  : 


_ 2xAT _ 

AX  AY  +  8X(2y  -  xy8Y)  ’ 
2yAX 

AX  AY  +  SY(2x  -  xySX )  ’ 


(20) 

(21) 


where  AX  =  (N  +  S)/2,  AY  =  (E  +  W)/2,  SX  =  (N-S)/2,  and  8Y  =  (E-W)/ 2.  These 
equations  are  nonlinear,  since  x  and  v  emerge  on  the  right-hand  side,  but  are  easily  iterated  to 
convergence.  Given  the  (x,y)  coordinates  of  the  midpoint,  we  relate  equation  (19)  to  the 
components  of  uat  the  cell  comers  to  estimate  the  velocity  at  the  midpoint.  The  midpoint 
velocity  is  changed  back  to  cell-corner  coordinates  using  the  inverse  of  equation  (17),  and  then 
the  corrected  velocity  in  equation  (18)  is  applied  to  find  the  departure  point.  With  this  correction, 
departure  points  for  a  velocity  field  varying  linearly  in  space  are  nearly  precise. 


5.2.2.13  Integrating  Fluxes 

In  the  next  step,  the  reconstructed  fields  are  integrated  over  the  departure  triangles  to  find  the 
total  fluxes  of  area,  volume,  and  energy  across  each  cell  edge.  Area  fluxes  are  simple  to  compute 
since  the  area  is  linear  in  x  and  y.  Given  a  triangle  with  vertices  xi=(xi,yi),  i  e  {1,2,3} ,  the 
triangle  area  is 


At=\ |(x2  “ Xj )( >3  - yl ) -  (y2  - Ti )(*3  (22) 

The  integral  /,  of  any  linear  function  / (r)  over  a  triangle  is  given  by 

1 1  =  ATf(x0),  (23) 
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where  x0  =  (x0 ,  v0 )  represents  the  triangle  midpoint, 


To  compute  the  area  flux,  the  area  is  evaluated  at  the  midpoint, 

a(x0)  =  ac+axx0+ayy0. 


(24) 


(25) 


and  multiplied  by  Ar .  By  convention,  northward  and  eastward  fluxes  are  positive,  while 
southward  and  westward  fluxes  are  negative. 

Equation  (23)  cannot  be  used  for  volume  fluxes,  since  the  reconstructed  volumes  are  quadratic 
functions  of  position.  (They  are  products  of  two  linear  functions,  area  and  thickness.)  The 
integral  of  a  quadratic  polynomial  over  a  triangle  necessitates  function  evaluations  at  three 
points, 


/2=4E/(xT  <26> 

3  i= 1 

where  x.  =  (x0  +x(.)/2  are  points  lying  halfway  between  the  midpoint  and  the  three  vertices.  The 

authors  of  [9]  use  this  fonnula  to  calculate  fluxes  of  the  product  pT ,  which  is  equivalent  to  ice 
volume.  Equation  (26)  does  not  work  for  ice  and  snow  energies,  which  are  cubic  functions- 
products  of  area,  thickness,  and  enthalpy.  Integrals  of  a  cubic  polynomial  over  a  triangle  can  be 
evaluated  using  a  four-point  formula  [42]: 


/3  At 


9_ 

16 


f(*  o)  + 


25 

48 


LaO 


(27) 


where  x.  =  (3x0  +2x.)/5  . 

To  evaluate  functions  at  specific  points,  products  of  the  form  a(x)h(x )  and  a(x)h(x)q(x)  must 
be  computed,  where  each  term  in  the  product  is  the  sum  of  a  cell-center  value  and  two 
displacement  terms.  This  calculation  may  be  sped  up  by  storing  some  sums  that  are  used 
repeatedly.  First,  weighted  ice  areas  are  computed  at  the  four  points  of  the  integral  formula  (27): 

a0  =  —^(ac+axx0+ayy0). 


ai  =  T^(ac  +  axxf  +  ayyi ),  /  e  {1,2,3}, 
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where  the  double  primes  have  been  dropped  from  the  jc.  .  Defined  next  is 

3 

Ga  = 

i= 0 
3 

i=0 

3 

%  = 

!=0 

which  is  used  to  compute  the  volume  fluxes: 


cr  ,  =  cr  h  +  cr  /z  +  cr  h  . 

ah  a  c  ax  x  ay  y 


Notice  that  oa,  crax,  and  <Jay  are  used  in  three  different  flux  computations:  ice  volume,  snow 
volume,  and  area- weighted  surface  temperature.  Defined  next  are 


cr 


cr 


cr 

ayy 

G axil 
^ ’ ayh 


’ 


i= 0 
3 


1=0 

3 


cr  h  +cr  h  +cr  h  , 

ax  c  axx  x  arp  y ? 

cr  h  +  cr  h  +  cr  h  . 

ay  c  axy  x  ayy  y 


The  sums  craxh  and  crayh  are  calculated  separately  for  ice  and  snow,  whereas  the  first  three  sums 

are  independent  of  the  material  being  transported.  Each  sum  is  utilized  repeatedly  if  there  are 
multiple  enthalpy  layers.  From  these  sums  the  energy  fluxes  are  computed: 

aahq  =  Gahilc  +  G axh^x  +(7ayh<ly ’ 

thus  finishing  the  flux  integrals  for  a  given  triangle.  To  compute  the  total  flux  across  a  cell  edge 
the  contributions  from  each  triangle  are  added. 
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5. 2. 2. 1. 4  Updating  State  Variables 

Finally,  these  fluxes  are  used  to  compute  new  values  of  the  state  variables  in  every  ice  category 
and  grid  cell.  The  new  fractional  ice  areas  a'in(i,j)  are  given  by 


’  /•  _  /•  • \  ,  FEa(j  l;j)  FEa(i,  j)  +  FNa(l,  j  1)  FNa(i,j")  /oq\ 

ai„Kl’J)-am\hJ)  + - 777^ -  UN 

Ahj) 

where  FEa(i,j )  and  Fm(i,j)  are  area  fluxes  across  the  east  and  north  edges,  respectively,  of  cell 
(z',y)  .  A(i,j )  is  the  grid  cell  area.  All  fluxes  added  to  one  cell  are  subtracted  from  a  neighboring 
cell;  thus  equation  (28)  conserves  total  ice  area. 

The  new  ice  volumes  and  energies  are  calculated  analogously.  New  thicknesses  are  provided  by 
the  ratio  of  volume  to  area,  and  enthalpies  by  the  ratio  of  energy  to  volume.  Tracer  monotonicity 
is  guaranteed  because 


h' 


ahdA 

J  A _ 

adA 

J  A 


where  h'  and  q'  are  the  new-time  thickness  and  enthalpy,  achieved  by  integrating  the  old-time 
ice  area,  volume,  and  energy  over  a  Lagrangian  departure  region  with  area  A.  In  other  words,  the 
new-time  thickness  and  enthalpy  are  weighted  averages  over  old-time  values,  with  non-negative 
weights  a  and  ah.  Therefore  the  new-time  values  must  lie  between  the  maximum  and  minimum 
of  the  old-time  values. 


5. 2. 2. 2  Transport  in  Thickness  Space 

Now  the  equation  for  ice  transport  in  thickness  space  due  to  thermodynamic  growth  and  melt  is 
solved  by, 

|4(/S)  =  0.  (29) 

ot  oh 

which  is  obtained  from  equation  (2)  by  ignoring  the  first  and  third  terms  on  the  right-hand  side. 
The  remapping  method  of  [25]  is  used,  in  which  thickness  categories  are  represented  as 
Lagrangian  grid  cells  with  boundaries  that  are  projected  ahead  in  time.  The  thickness  distribution 
function  g  is  estimated  as  a  linear  function  of  h  in  each  displaced  category  and  is  then  remapped 
onto  the  initial  thickness  categories.  This  method  is  numerically  smooth  (unlike  schemes  that 
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treat  g{h )  like  a  set  of  delta  functions)  and  is  not  overly  diffusive.  It  can  be  seen  as  a  ID 
simplification  of  the  2D  incremental  remapping  scheme  described  above. 

The  first  computation  is  the  displacement  of  category  boundaries  in  thickness  space.  Assume  that 
at  time  m  the  ice  areas  a”  and  mean  ice  thicknesses  /?"'  are  known  for  each  thickness  category. 
(For  now  the  subscript  i  that  distinguishes  ice  from  snow  must  be  omitted.)  A  thennodynamic 
model  (Section  5. 2.2. 5)  is  used  to  calculate  the  new  mean  thicknesses  h’’,+1  at  time  m  +1.  The 

time  step  must  be  small  enough  that  trajectories  do  not  cross;  i.e.,  h'"+l  <  /zjjj1  for  each  pair  of 
adjacent  categories.  The  growth  rate  at  h  =  hn  is  given  by  fn  =  (/?"'"  —  h”)/At .  By  linear 
interpolation  growth  rate  F„  is  estimated  at  the  upper  category  boundary  //„: 

F„=f,+TBlZTL(H,-h,). 

hn+l~hn 

If  an  or  an+l  =  0  ,  Fn  is  set  to  the  growth  rate  in  the  nonzero  category,  and  if  a„  =  a„+ i  =0,  we  set 
Fn  =0.  The  temporary  displaced  boundaries  for  n  =1  to  N  -  1  are  given  by 

H*  =  H+FAt. 

n  n  n 

The  boundaries  cannot  be  displaced  by  more  than  one  category  to  the  left  or  right;  in  other 
words,  Hn]  <  H*n  <  Hn+1  is  required.  Without  this  requirement  there  would  need  to  be  a  general 
remapping  rather  than  an  incremental  remapping,  at  the  cost  of  added  complexity. 

Next  g(h)  is  constructed  in  the  displaced  thickness  categories.  The  ice  areas  in  the  displaced 
categories  area”+1  =  <?"' ,  because  area  is  conserved  following  the  motion  in  thickness  space  (i.e., 
during  vertical  ice  growth  or  melting).  The  new  ice  volumes  arc  v’” 1 1  =  (anhn)m+1  =  .  For 

brevity,  define  //,  =  H*n  ,  and  II R  =  H*n  and  drop  the  time  index  m  +1.  A  continuous  function 
g(h)  is  built  in  to  each  category  so  that  the  total  area  and  volume  at  time  m  +1  are  a„  and  v„, 
respectively: 


(30) 


r%r<a  =  V,.  (31) 

jhl 

The  simplest  polynomial  that  can  satisfy  both  equations  is  a  line.  It  is  easy  to  change  coordinates, 
writing g(r/)  =  g,//  +  g0 ,  where  r/  =  h-HL  and  the  coefficients  g0  and  g,  are  to  be  determined. 
Then  equations  (30)  and  (3 1)  may  be  written  as 
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8i  y  +  go^  =an, 
rjl  rj\ 


where  tjr  =  HR-  II ,  and  t]n  =  hn  -  II , .  The  solution  to  these  equations  is 


^0 


_  6«„ 


vl 


2?Ir 


n„ 


s  i  =  — 

Vr 


Xlan[n 

In  rs 

V  z 


Since  g  is  linear,  its  minimum  and  maximum  values  lie  at  the  boundaries,  rj  =  0  and///( 


(32) 

(33) 


g(0)  = 


6  aj  2  rjR 


Vr  \  3 


>1, 


8o> 


g(7*)  =  % 
rR 


}In 


Vr 


(34) 

(35) 


Equation  (34)  implies  that  g(0)  <  0  when  r/n  >  2jjr/3  ,  i.e.,  when  hn  sits  in  the  right  third  of  the 
thickness  range  ( Hl,Hr ).  Similarly,  equation  (35)  means  that  g(rjR)<  0  when  tjn  <  V3 ,  i-e., 
when  hn  is  in  the  left  third  of  the  range.  Because  negative  values  of  g  are  unphysical,  a  different 
solution  is  used  when  h„  lies  outside  the  central  third  of  the  thickness  range.  If  hn  is  in  the  left 
third  of  the  range,  a  cutoff  thickness  is  identified  as  Hc  =3 hn  -  2 HL,  and  sets  g  =0  between  1 1(  and 
Hr.  Equations  (32)  and  (33)  are  then  valid  with  r/R  redefined  as  He  -  HL.  And  if  hn  is  in  the  right 
third  of  the  range,  define  He  =3 hn  -  2 Hr  and  set  g  =0  between  ///  and  He-  In  this  case,  equations 
(32)  and  (33)  apply  with  t)R  =  HR-HC  and  //„  =  hn  -  H c . 

Figure  4  demonstrates  the  linear  reconstruction  of  g  for  the  simple  cases  HL  =0,  Hr  =1,  a„  =1, 
and  hn  =  0.2,  0.4,  0.6,  and  0.8.  Note  that  g  slopes  downward  (gj  <  0)  when  hn  is  less  than  the 
midpoint  thickness,  ( HL  +  Hr)/ 2=1/2,  and  upward  when  h„  surpasses  the  midpoint  thickness.  For 
h„  =0.2  and  0.8,  g  =0  over  part  of  the  range. 
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Figure  4:  Linear  approximation  of  the  thickness  distribution  function  g(h)  for  an  ice 
category  with  left  boundary  HL  =0,  right  boundary  HR  =1,  fractional  area  a„  =1,  and  mean 
ice  thickness  hn  =  0.2,  0.4,  0.6,  and  0.8. 

Finally,  the  thickness  distribution  is  remapped  to  the  original  boundaries  by  transferring  area  and 
volume  between  categories.  The  ice  area  A a„  and  volume  Av„  are  calculated  between  each 
original  boundary  Hn  and  displaced  boundary//* .  If,  H*n  >  Hn ,  ice  changes  from  category  n  to  n 
+ 1 .  The  area  and  volume  transferred  are 

Aa"  =  \l!"gdh'  (36) 

nn 

Av„  =  hgdh .  (37) 

n  n 

if  K  <  Hn,  ice  area  and  volume  are  shifted  from  category  n  +1  to  n  using  equations  (36)  and 
(37)  with  the  limits  of  integration  reversed.  To  evaluate  the  integrals,  coordinates  are  changed 
from  h  to  rj  =  h-  Hi,  where  Hl  is  the  left  limit  of  the  range  over  which  g  >  0,  and  write  g(rj)  using 
equations  (32)  and  (33).  This  is  how  the  new  areas  a„  and  volumes  v„  are  acquired  between  the 
original  boundaries  //„_;  and  Hn  in  each  category.  The  new  thicknesses,  hn  =  vjan ,  are 

guaranteed  to  be  within  the  range  (Hn-i,  Hn).  If  g  =0  in  the  area  of  a  category  that  is  remapped  to 
a  neighboring  category,  no  ice  is  transferred. 

Other  conserved  quantities  are  transported  in  proportion  to  the  ice  volume  Avin.  (Now  using  the 
subscripts  i  and  s  to  distinguish  ice  from  snow.)  For  example,  the  transferred  snow  volume  is 
Avs.„  =  vsn(Avm/vin),  and  the  transferred  ice  energy  in  layer  k  is  A <?,„/,  =  eink(Avin/vin). 


n-  0.2 

h-  0.4 
h-  0.6 
h-  0.8 
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The  left  and  right  boundaries  of  the  domain  necessitate  special  treatment.  If  ice  is  accumulating 
in  open  water  at  a  rate  F(h  the  left  boundary  Ho  is  shifted  right  by  FoAt  before  g  is  constructed  in 
category  1,  then  reset  to  zero  after  the  remapping  is  finished.  Then  new  ice  is  added  to  the  grid 
cell,  conserving  area,  volume,  and  energy.  If  ice  cannot  grow  in  open  water  (if  the  sea  water  is 
too  warm  or  the  net  surface  energy  flux  is  downward),  H0  is  fixed  at  zero,  and  the  growth  rate  at 
the  left  boundary  is  approximated  as  Fq  =fi.  If  F0  <  0,  the  area  of  ice  thinner  than  A  ho  =  -FoAt  is 
added  to  the  open  water  area  a0,  allowing  the  ice  and  snow  volume  and  energy  to  remain 
unchanged.  The  area  of  new  open  water  is 


f  A/i0 

Aao=J()  gdh. 

The  right  boundary  HN  is  not  fixed  but  varies  with  hu,  the  mean  ice  thickness  in  the  thickest 
category.  Given  /z,v,  set  HN  =3 hN  -  2//V-/,  which  guarantees  that  g(h)  >  0  for  //.v_ /  <  h  <  HN  and 
g(h)= 0  for  h  >  IIh< .  No  ice  crosses  the  right  boundary. 

If  the  ice  growth  or  melt  rates  in  a  given  grid  cell  are  too  big,  the  thickness  remapping  scheme 
will  not  function.  Instead,  the  thickness  categories  in  that  grid  cell  are  treated  as  delta  functions 
following  [6],  and  categories  outside  their  set  boundaries  are  merged  with  neighboring  categories 
as  needed.  For  time  steps  of  less  than  a  day  and  category  thickness  ranges  of  10  cm  or  more,  this 
simplification  is  rarely  needed,  if  ever. 


5.2. 2. 3  Mechanical  Redistribution 


The  last  term  on  the  right-hand  side  of  equation  (2)  is  i//,  which  defines  the  redistribution  of  ice  in 
thickness  space  due  to  ridging  and  other  mechanical  processes.  The  mechanical  redistribution 
scheme  in  PIPS  3.0  is  based  on  [43],  [35],  [16],  [27]  and  [13].  This  scheme  transforms  thinner 
ice  to  thicker  ice  and  is  applied  after  horizontal  transport.  When  the  ice  is  converging,  enough  ice 
ridges  to  guarantee  that  the  ice  area  does  not  surpass  the  grid  cell  area. 


First  specify  the  participation  function:  the  thickness  distribution  ap(h)  =  b(h)g(h)  of  the  ice 

participating  in  ridging.  (“Ridging”  is  used  here  as  shorthand  for  all  modes  of  mechanical 
redistribution.)  The  weighting  function  b(h)  prefers  ridging  of  thin  ice  and  closing  of  open  water 
to  ridging  of  thicker  ice.  There  are  two  options  for  the  form  of b(h) .  If  krdgpartic  =  0  in  the 

namelist,  then  following  [43],  we  set 


b(h)  = 


2 

¥ 

o 


(i 


G(A) 

G 


ifG(h)  <  G* 
otherwise 


(38) 
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where  G(h )  is  the  fractional  area  covered  by  ice  thinner  than  h,  and  G*  is  an  empirical  constant. 
Integrating  ap(h)  between  category  boundaries  //„_/  and  Hn,  the  mean  value  of  ap  in  category  n 
is  obtained  by 


*Pn 


=  ^(G„-Gnl) 


1- 


G^+G. 


2G 


(39) 


where  aPn  is  the  ratio  of  the  ice  area  ridging  (or  open  water  area  closing)  in  category  n  to  the 
total  area  ridging  and  closing,  and  Gn  is  the  total  fractional  ice  area  in  categories  0  to  N.  Equation 
(39)  applies  to  categories  with  Gn  <  G*.  If  G„.i  <  G*  <Gn,  equation  (39)  is  valid  with  G* 
replacing  Gn,  and  if  Gn.i  >  G*,  then  aPn  =0.  If  the  open  water  fraction  ao  >G*,  no  ice  will  ridge, 

because  “ridging”  simply  reduces  the  area  of  open  water.  As  in  [43],  G*  is  set  to  =0. 15. 

If  the  spatial  resolution  is  too  fine  for  a  given  time  step  At,  the  weighting  function  (38)  can 
support  numerical  instability.  For  At  =  1  hour,  resolutions  finer  than  Ax  ~  1 0  km  are  usually 
unstable.  The  instability  comes  from  feedback  between  the  ridging  scheme  and  the  dynamics 
through  the  ice  strength.  If  the  strength  changes  considerably  on  time  scales  less  than  At,  the 
EVP  solution  becomes  inaccurate  and  sometimes  oscillatory.  Consequently,  the  fields  of  ice 
area,  thickness,  velocity,  strength,  divergence,  and  shear  can  become  noisy  and  unphysical. 

A  more  stable  weighting  function  was  suggested  by  [27]: 


=  exp[-G(/i)/a*] 

1  '  fl*[l-exp(-l/fl*)] 


When  integrated  between  category  boundaries,  equation  (40)  implies 

=  exp(-G„/a* )  -  exp(-G„  Ja  ) 
l-exp(-l/<7*) 


(40) 


(41) 


This  weighting  function  is  applied  if  krdgpartlc  =  1  in  the  namelist.  From  equation  (40),  the 

mean  value  of  G  for  ice  involved  in  ridging  is  a*,  as  compared  to  G*  1 3  for  equation  (38).  For 
standard  ice  thickness  distributions,  setting  a*  =  0.05  with  krdgpartic  =  1  generates 

participation  fractions  similar  to  those  given  by  G*  =  0.15  withkrdgpartic  =  0  .  See  [27]  for  a 
detailed  comparison  of  these  two  participation  functions. 

Thin  ice  is  converted  to  thick  ridged  ice  in  a  way  that  diminishes  the  total  ice  area  while 
conserving  ice  volume  and  energy.  There  are  two  namelist  options  for  redistributing  ice  among 
thickness  categories.  Ifkrdgredist  =  0 ,  ridging  ice  of  thickness  h„  creates  ridges  whose  area  is 
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distributed  uniformly  between  Hmin  =  2 hn  and  Hmdx  =  2 y]H*hn  ,  as  in  [16].  The  default  value  of 

H*  is  25  m.  Observations  suggest  that  II  =  50  m  gives  a  better  fit  to  first-year  ridges  [4], 
although  the  lower  value  may  be  appropriate  for  multiyear  ridges  [13].  The  ratio  of  the  mean 
ridge  thickness  to  the  thickness  of  ridging  ice  is  kn  =  (//min  +  f/max)/(2/z;() .  If  the  area  of  category 

n  is  reduced  by  ridging  at  the  rate  r„,  the  area  of  thicker  categories  matures  simultaneously  at  the 
rate  rn/kn.  Thus  the  net  rate  of  area  loss  due  to  ridging  of  ice  in  category  n  is  rn{  1  -  1/A:,,).  The 
ridged  ice  area  and  volume  are  allocated  amid  categories  in  the  thickness  range (Hmin,Hmax)  •  The 
fraction  of  the  new  ridge  area  going  to  category  m  is 


L 


area  Hr  ^ L 


TT  _  TT 

1  max  min 


(42) 


where  HL  =  max(//n  )  and  II R  =  min(//m,  //max ) .  The  fraction  of  the  ridge  volume  going 

to  category  m  is 


(Hr)2  ~{H l)2 
(Hmj2-(Hmmr 


(43) 


This  uniform  redistribution  function  has  a  tendency  to  not  produce  enough  ice  in  the  3—5  m 
range  and  too  much  ice  thicker  than  10  m  [4].  Observational  data  show  that  the  ITD  of  ridges  is 
better  approximated  by  a  negative  exponential.  Setting  krdgredist  =  1  gives  ridges  with  an 
exponential  ITD  [27]: 


gR(h)KQxp[-(h-Hmm)U]  (44) 

for  h  >=  Hmm  ,  with  gR  ( h )  =  0  for  h  <  Hmm .  Here,  the  variable  A  is  an  empirical  e-folding  scale 
and  Hmm  =  2 hn  (where  hn  represents  the  thickness  of  ridging  ice).  It  is  assumed  that  A  =  /uh)'2 , 
where  //  is  a  tunable  parameter  with  units  of  m1 2 .  Therefore  the  mean  ridge  thickness  increases 
in  proportion  to  hlJ2 ,  as  in  [16].  The  default  value  //  =  4.0  m 1/2  assigns  A  in  the  range  1—4  m  for 
most  ridged  ice.  Ice  strengths  with  ju  =  4.0  m1 2  and  krdgredist  =  1  are  nearly  comparable  to  the 
strengths  with  H*  =  50  and  krdgredist  =  0  . 

From  equation  (44)  it  can  be  shown  that  the  fractional  area  going  to  category  m  as  a  result  of 
ridging  is 


fr  =  exp[-(/A„,_1  -  Hmin)/A\ -  exp[-(//m -Hmm)/A\.  (45) 


The  fractional  volume  available  to  category  m  is 
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(Hm-i  +^)cxp[-(//„_, -//mm )/i]-(//„+i)exp[-(//„  -Hmin)/JL\ 


(46) 


Equations  (45)  and  (46)  replace  equations  (42)  and  (43)  whenkrdgredlst  =  1 . 

Internal  ice  energy  is  moved  between  categories  in  proportion  to  ice  volume.  Snow  volume  and 
energy  are  moved  in  the  same  way,  except  that  a  fraction  of  the  snow  could  be  deposited  in  the 
ocean  instead  of  added  to  the  new  ridge. 

The  net  area  removed  through  ridging  and  closing  is  a  function  of  the  strain  rates.  Let  Rna  equal 
the  net  rate  of  area  loss  for  the  ice  pack  (i.e.,  the  rate  of  open  water  area  closing,  in  addition  to 
the  net  rate  of  ice  area  loss  due  to  ridging).  Following  [13],  Rncl  is  given  by 

^net  =  (A- 1  Dd  I)  -  min(Z)Z) ,  0),  (47) 


where  Cs  is  the  fraction  of  shear  dissipation  energy  that  contributes  to  ridge-building,  Dn 
represents  the  divergence,  and  A  is  a  function  of  the  divergence  and  shear.  These  strain  rates  are 
calculated  by  the  dynamics  scheme.  The  default  value  of  Cs  is  0.25. 

Define  Rtot  =  ^  78 .  This  rate  is  related  to  Rna  by 


f? 


apo 


Pn 


(1-1^) 


R, 


(48) 


Given  Rnet  from  equation  (47),  equation  (48)  is  used  to  compute  RUtl .  Next  the  area  ridged  in 
category  n  is  shown  by  a  rn  =  rAt ,  where  r  =  aPnRtot .  The  area  of  new  ridges  is  ajkn ,  and  arnhn 

represents  the  volume  of  new  ridges  (because  volume  is  conserved  during  ridging).  The  ridging 
ice  from  category  n  is  removed  and  we  use  equations  (42)  and  (43)  or  equations  (45)  and  (46)  to 
redistribute  the  ice  amid  thicker  categories. 

In  some  instances  the  ridging  rate  in  a  given  thickness  category  n  may  be  substantial  enough  to 
ridge  the  total  area  in  that  category  during  a  time  interval  less  than  At.  In  this  instance  Rtot  is 
reduced  to  the  value  that  accurately  ridges  an  area  a„  during  At.  Following  each  ridging  iteration, 
the  total  fractional  ice  area  a,  is  computed.  If  at  >  1,  the  ridging  is  repeated  with  a  value  of  Rnet 
sufficient  to  yield  a,  =1 . 

The  ice  strength  P  may  be  computed  in  either  of  two  ways.  If  the  namelist  parameter 
krdgredlst  =  0 ,  the  strength  is  given  by  [15]: 
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P  =  P*h  exp[-C ( 1  -  a,. )] , 


(49) 


where  P*  =  27,  500  N/m  and  C  =  20  are  empirical  constants  and  h  is  the  mean  ice  thickness. 
Alternatively,  setting  krdgredist  =  1  provides  an  ice  strength  closely  related  to  the  ridging 

scheme.  Following  [35],  the  strength  is  known  to  be  proportional  to  the  change  in  ice  potential 
energy  A Ep  per  unit  area  of  compressive  deformation.  Given  the  uniform  ridge  ITDs 
( krdgredlst  =  0 ),  notice  that 


p=crcrj£ 


.a  h2+^~ 

UPnUn  ^  r 

k, 


n  V 


(Err -wr) 

T/rrm  ax  rrminx 

3(H„  ~En  ) 


3  A 


(50) 


where  CP  =  (gl2)(pj pw)(pw  -pt)  ,  =  Rtot/Rnet>  1  from  equation  (48),  and  Cf  represents  an 

empirical  parameter  that  accounts  for  frictional  energy  dissipation.  Following  [13],  we  set  C/  = 
17.  The  primary  tenn  in  the  summation  is  the  potential  energy  of  ridging  ice,  and  the  second, 
greater  tenn  is  the  potential  energy  of  the  resulting  ridges.  The  factor  of  p  is  included  since  ap„  is 
nonnalized  with  respect  to  the  total  area  of  ice  ridging,  not  the  net  area  removed.  Remember  that 
more  than  one  unit  area  of  ice  must  be  ridged  to  reduce  the  net  ice  area  by  one  unit.  For 
exponential  ridge  ITDs  (krdgredlst  =  1),  the  ridge  potential  energy  is  adapted: 


n= 1 


-«n,*.i+TL(Ai„+2ff„i,3  +  2a2) 


(51) 


The  energy-based  ice  strength  given  by  equations  (50)  or  (51)  is  more  physically  realistic 
than  the  strength  produced  by  equation  (49).  But  the  use  of  equation  (49)  is  less  likely  to  pennit 
numerical  instability  at  a  given  resolution  and  time  step.  See  [27]  for  more  details. 


5. 2. 2. 4  Dynamics 

The  elastic-viscous-plastic  (EVP)  model  represents  an  alteration  of  the  standard  viscous-plastic 
(VP)  model  for  sea  ice  dynamics  by  [15].  It  reduces  to  the  VP  model  at  temporal  scales 
associated  with  the  wind  forcing,  while  at  shorter  time  scales  the  modification  process  takes 
place  by  a  numerically  more  efficient  elastic  wave  mechanism.  While  retaining  the  basic  physics, 
this  elastic  wave  version  leads  to  a  fully  explicit  numerical  scheme  that  greatly  enhances  the 
model’s  computational  efficiency. 

The  EVP  sea  ice  dynamics  model  is  well  documented  in  [18],  [17],  [19]  and  [20],  Simulation 
results  and  performance  of  the  EVP  model  have  been  evaluated  against  the  VP  model  in  realistic 
simulations  of  the  Arctic  [21].  The  equations  are  summarized  below  but  it  is  recommended  that 
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the  reader  refer  back  to  the  above  references  for  details.  The  numerical  implementation  in  this 
code  release  is  that  of  [19],  [20]. 

The  force  balance  per  unit  area  in  the  ice  pack  is  set  by  a  two-dimensional  momentum  equation 
[15],  obtained  by  integrating  the  3D  equation  through  the  thickness  of  the  ice  in  the  vertical 
direction: 


du  _ 

m  —  =  V  •  a  +  r  +  r 
dt 


k  x  mfu  -  mgVHo , 


(52) 


where  m  is  the  combined  mass  of  ice  and  snow  per  unit  area  and  fa  and  fw  represent  wind  and 
ocean  stresses,  respectively.  The  strength  of  the  ice  is  given  by  the  internal  stress  tensor  cr  ,  and 

the  other  two  terms  on  the  right  side  are  stresses  due  to  Coriolis  effects  and  the  sea  surface  slope. 
The  parameterization  for  the  wind  and  ice-ocean  stress  tenns  must  have  the  ice  concentration  as 
a  multiplicative  factor  to  be  consistent  with  the  formal  theory  of  free  drift  in  low  ice 
concentration  areas.  A  thorough  explanation  of  the  issue  and  its  continuum  solution  is  given  in 
[20]  and  [8], 

For  simplicity  the  stress  tensor  a  is  formulated  in  terms  of  crl  =  an  +  a22,  a2  =  an  -a22  and  the 
divergence,  Dd  The  stress  tensor  is  also  formulated  in  tenns  of  the  horizontal  tension  and 
shearing  strain  rates,  Dt  and  Ds  respectively.  The  internal  stress  tensor  is  determined  from  a 
standardized  version  of  the  VP  constitutive  law, 


1  da,  a ,  P 
- -  H - —  + -  = 

E  dt  2C,  2C, 

±Al+^l  =  Drt 

E  dt  2?7 
J_da11+a11  =  ]_D 

E  dt  2rj  2 


(53) 

(54) 

(55) 


where 
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dt 

Ds 


c 

v 


£n  +£-22’ 
^11  —£22’ 


2s, 


12’ 


GU  dll  ■ 

— -  + — - 
dxf  fix,.  . 

V  1 J 


_P_ 

2A  ’ 
P 

2A^ 


(56) 

(57) 

(58) 


A 


D2D+\e2(D2T 


and  P  is  a  function  of  the  ice  thickness  and  concentration,  explained  in  Section  5. 2.2. 3.  The 
dynamics  component  utilizes  a  “replacement  pressure”  (see  [14],  for  example),  which  functions 
to  prevent  residual  ice  motion  due  to  spatial  variations  of  P  when  the  rates  of  strain  are  exactly 
zero. 


Many  adjustments  have  been  made  to  the  EVP  model  since  its  original  release.  In  the  previous 
version,  the  viscosities  were  held  fixed  while  the  stress  and  momentum  equations  were  subcycled 
with  the  smaller  time  step  dte.  The  reason  for  applying  the  EVP  model  in  this  way  was  to 
reproduce  the  results  of  the  original  VP  model  as  closely  as  possible.  When  solved  with  time 
steps  of  several  hours  or  more,  the  VP  model  goes  through  a  linearization  error  associated  with 
the  viscosities,  which  are  lagged  over  the  time  step  [17].  This  led  to  principal  stress  states  that 
were  widely  spread  outside  the  elliptical  yield  curve  in  both  models  [21].  This  problem  has  been 
addressed  by  updating  the  viscosities  during  the  subcycling,  so  that  the  complete  dynamics 
component  is  subcycled  within  the  time  step.  Taken  alone,  this  modification  would  require  an 
increased  number  of  operations  to  compute  the  viscosities. 

Nonetheless,  the  new  dynamics  component  is  roughly  as  efficient  as  the  earlier  version  due  to  a 
change  in  the  definition  of  the  elastic  parameter  E.  E  is  now  described  in  terms  of  a  damping 
timescale  T  for  elastic  waves,  A te<T  <  At ,  as 


E  = 


T’ 


where  T  =  EoAt  and  Eo  (eye)  is  a  tunable  parameter  less  than  one,  as  before.  The  stress 
equations  (53-55)  become 
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da ,  a,  P 
- -  H - —  + - 

8t  2  T  2  T 
da2  e2a2 
dt  2  T 
dan  |  e  al2 
dt  2  T 

All  coefficients  on  the  left  side  are  constant  with  the  exception  of  P,  which  changes  only  on  the 
longer  time  step  At.  This  alteration  compensates  for  the  diminished  efficiency  of  including  the 
viscosity  terms  in  the  subcycling.  (Note  that  the  viscosities  do  not  appear  explicitly.)  Choices  of 
the  parameters  used  to  describe  E,  T and  A te  are  discussed  in  the  PIPS  3.0  User’s  Manual  [2]. 

A  different  discretization  of  the  stress  terms  da^/dx.  in  the  momentum  equation  is  now  used, 

which  allows  the  discrete  equations  to  be  derived  from  the  continuous  equations  written  in 
curvilinear  coordinates.  In  this  way,  metric  terms  associated  with  the  curvature  of  the  grid  were 
integrated  into  the  discretization  explicitly.  No  longer  in  use  is  the  “triangle  discretization,”  in 
which  the  strain  rates  and  stresses  were  constant  over  each  of  four  subtriangles  in  each  grid  cell, 
and  instead  a  bilinear  approximation  for  the  velocities  and  stresses  is  used.  Details  pertaining  to 
the  spatial  discretization  are  given  in  [19]. 

The  momentum  equation  is  discretized  in  time  as  follows.  First,  to  clarify,  the  two  components 
of  equation  (52)  are 


2T  A 


-D 


D  5 


2TA 

P 


AT  A 


Dt, 


Ds. 


mJt  =  d~^~  +Tax+aicwPw\uw-u\[(Uw-u)cos0-(Vw-v)sm0]  +  mfv-mg ^ 

j  j 

dv  -  +  Tay  +  a  fKP,  pw  - 11 1  [( uw  -  U  )  sin  61  -  ( Vw  -  v)  cos  6>]  -  mfu  -  mg 


m 


dt 


In  the  code,  vrel  =  aicwp\Uw-u  \,  where  k  denotes  the  subcycling  step.  The  following 
equations  show  the  time  discretization  and  illustrate  some  of  the  other  variables  used  in  the  code. 


Jfl 

—  +  vrel  cos  6 
At  j 


^£j^+i  P)TT 

u A+1  -  ( mf  +  vrel  sin  0)  vk+l  =  — —  +  rax  -  mg  — 1  +  vrel  ( U  cos  0  -  V  sin  0)  +  —  uk , 

x,  dx  - -  At 


ccb 


J 

strintx 


forcex 


ccb 


[mf  +  vrel  sin  !l  +[  -^-  +  vrelcos0 


'ATT 

vm  =  — —  +  r  -  mg  - — 1  +  vrel  ( Uw  sin  6  +  V  cos  6 )  +  —  v* , an(^ 
dx  ~ - - •  At 


forcey 


watery 
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vrel  -waterx  (y)  =  taux  (y) .  This  system  of  equations  is  solved  analytically  for  uk+1  and 
vA+1 .  When  the  subcycling  is  complete  for  each  (thermodynamic)  time  step,  the  ice-ocean  stress 
must  be  constructed  from  taux(y)  and  the  terms  containing  vrel  on  the  left  side  of  the 
equations.  This  is  done  in  subroutine  ev p_finish. 


5. 2. 2. 5  Thermodynamics 

The  thennodynamic  sea  ice  model  is  based  on  [31]  and  [7],  and  is  described  more  fully  in  [24]. 
For  each  thickness  category  the  model  calculates  changes  in  the  ice  and  snow  thickness  and 
vertical  temperature  profile  resulting  from  radiative,  turbulent,  and  conductive  heat  fluxes.  The 
ice  has  a  temperature-dependent  specific  heat  to  reproduce  the  effect  of  brine  pocket  melting  and 
freezing. 

Each  thickness  category  n  in  each  grid  cell  is  handled  as  a  horizontally  uniform  column  with  ice 
thickness  hm  =  v,n/aw  and  snow  thickness  hsn  =  vsn/ain.  (Henceforth  the  category  index  n  is 
omitted.)  Each  column  is  divided  into  /V,  ice  layers  of  thickness  Ah,  =  h,/N,  and,  if  snow  is 
present,  a  single  snow  layer.  (Allowing  for  multiple  snow  layers  is  possible  in  future  versions  of 
PIPS  3.0.)  The  surface  temperature  (i.e.,  the  temperature  of  ice  or  snow  at  the  interface  with  the 
atmosphere)  is  Tsfi  which  cannot  exceed  0°C.  The  temperature  at  the  midpoint  of  the  snow  layer 
is  Ts,  and  the  midpoint  ice  layer  temperatures  are  7)*,  where  k  ranges  from  1  to  Nu  The 
temperature  at  the  ice  bottom  is  held  at  Tf,  the  freezing  temperature  of  the  ocean  mixed  layer.  All 
temperatures  are  in  degrees  Celsius  unless  otherwise  noted. 

The  vertical  salinity  profile  is  set  and  is  unchanging  in  time.  The  snow  is  assumed  to  be  fresh, 
and  the  midpoint  salinity  in  each  ice  layer  is  given  by 

(59) 

where  z  =  (k- 1/2 )/Ni ,  Smax  =3.2  psu.  The  variables  a  =0.407  and  b  =0. 573  are  determined  from 

a  least-squares  fit  to  the  salinity  profile  seen  in  multiyear  sea  ice  by  [36].  This  profile  varies  from 
S  =0  at  the  top  surface  (z  =0)  to  5’  =  Smax  at  the  bottom  surface  (z  =1)  and  is  similar  to  that  used 
by  [31].  Equation  (59)  is  quite  accurate  for  ice  that  has  drained  at  the  top  surface  due  to  summer 
melting.  It  is  not  a  good  approximation  for  cold  first-year  ice,  though,  which  has  a  more 
vertically  uniform  salinity  since  it  has  not  yet  drained.  However,  the  effects  of  salinity  on  heat 
capacity  are  negligible  for  temperatures  well  below  freezing,  so  the  salinity  error  does  not  lead  to 
significant  temperature  errors. 

Each  ice  layer  has  an  enthalpy  q defined  as  the  negative  of  the  energy  necessary  to  melt  a  unit 
volume  of  ice  and  raise  its  temperature  to  0°C.  Due  to  internal  melting  and  freezing  in  brine 
pockets,  the  enthalpy  of  the  ice  depends  on  the  brine  pocket  volume  and  is  a  function  of 
temperature  and  salinity.  Since  the  salinity  is  prescribed,  there  is  a  one-to-one  relationship 
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between  temperature  and  enthalpy.  Snow  enthalpy  qs  is  also  defined,  which  depends  on 
temperature  alone.  Enthalpy  equations  are  derived  in  Section  5.2. 2. 5. 3. 

Given  surface  forcing  at  the  atmosphere-ice  and  ice-ocean  interfaces  in  addition  to  the  ice  and 
snow  thicknesses  and  temperatures/enthalpies  at  time  m,  the  thermodynamic  model  advances 
these  quantities  to  time  m  +1.  The  computation  proceeds  in  two  steps.  First  to  be  solved  is  a  set 
of  equations  for  the  new  temperatures,  as  discussed  in  Section  5. 2. 2. 5. 2.  Next  is  the  calculation 
of  the  melting,  if  any,  of  ice  or  snow  at  the  top  surface,  and  the  growth  or  melting  of  ice  at  the 
bottom  surface,  as  described  in  Section  5. 2. 2. 5. 3.  We  begin  by  defining  the  surface  forcing 
parameterizations,  which  are  closely  related  to  the  ice  and  snow  surface  temperatures. 


5. 2.2.5. 1  Thermodynamic  Surface  Forcing 

The  net  energy  flux  from  atmosphere  to  ice  (with  all  fluxes  defined  as  positive  downward)  is 


F0  =  Fs  +F. +FLl+FLT+V-a)(l-i0)Fs 


where  Fs  is  the  sensible  heat  flux,  F\  is  the  latent  heat  flux,  FL^  is  the  incoming  longwave  flux,  F^ 
is  the  outgoing  longwave  flux,  Fsw  is  the  incoming  shortwave  flux,  a  represents  the  shortwave 
albedo,  and  io  is  the  fraction  of  absorbed  shortwave  flux  that  penetrates  the  ice. 

The  albedo  depends  on  the  thickness  and  temperature  of  ice  and  snow  and  on  the  spectral 
distribution  of  the  incoming  solar  radiation.  Albedo  parameters  have  been  selected  to  fit 
observations  from  the  Surface  Heat  Budget  of  the  Arctic  Ocean  (SHEBA)  field  experiment.  For 
Tsf  <  -1°C  and  ht  >  0.5m,  the  bare  ice  albedo  is  0.78  for  visible  wavelengths  (<  700  nm)  and 
0.36  for  near  IR  wavelengths  (>  700  nm).  As  ht  decreases  from  0.5  m  to  zero,  the  ice  albedo 
decreases  efficiently  (through  use  of  an  arctangent  function)  to  the  ocean  albedo,  0.06.  The  ice 
albedo  in  both  spectral  bands  drops  by  0.075  as  T,/ rises  from  -1°C  to  0°C.  The  albedo  of  cold 
snow  (Tsf<  -1°C)  is  0.98  for  visible  wavelengths  and  0.70  for  near  IR  wavelengths.  The  visible 
snow  albedo  falls  by  0.10  and  the  near  IR  albedo  by  0.15  as  7)/- increases  from  -1°C  to  0°C.  The 
total  albedo  is  an  area-weighted  average  of  the  ice  and  snow  albedos,  where  the  fractional  snow- 
covered  area  is 


f 

J  sn 


h  +  h 


snowpatch 


and  hSnowpatch  =0.02  m.  The  envelope  of  albedo  values  is  shown  in  Figure  5. 
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Figure  5:  Albedo  as  a  function  of  ice  thickness  and  temperature  for  the  two  extrema  in 
snow  depth.  Maximum  snow  depth  is  calculated  based  on  Archimedes’  Principle  for  the 
given  ice  thickness.  These  curves  symbolize  the  envelope  of  potential  albedo  values. 

The  net  absorbed  shortwave  flux  is  Fswabs  =  la  -0Cj)Fswir.,  where  the  summation  spans  four 

radiative  groups  (direct  and  diffuse  visible,  direct  and  diffuse  near  IR).  The  flux  penetrating  the 
ice  is  given  by/0  =  i0Fswabs ,  where  io  =0.70  (1  -  fsnow)  for  visible  radiation  and  i0  =0  for  near  IR. 
Radiation  penetrating  the  ice  is  attenuated  according  to  Beer’s  Law: 

I(z)  =  I0  exp (-/c,.z),  (60) 

where  I(z)  is  the  shortwave  flux  that  extends  to  depth  z  beneath  the  surface  without  being 

i 

absorbed,  and  Kt  is  the  bulk  extinction  coefficient  for  solar  radiation  in  ice,  set  at  1.4  in'  for 
visible  wavelengths  [12].  A  fraction  exp  (-/c,7i/)  of  the  penetrating  solar  radiation  passes  through 
the  ice  to  the  ocean  ( F^).  Although  incoming  shortwave  and  longwave  radiation  are  received 

from  the  atmosphere,  outgoing  long-wave  radiation  and  the  turbulent  heat  fluxes  are  derived 
quantities.  Outgoing  longwave  assumes  the  standard  blackbody  form ,Fl^  =  £ct{t^  ,  where 
s  =  0.95  is  the  emissivity  of  snow  or  ice,  cr  is  the  Stefan-Boltzmann  constant  and  is  the 
surface  temperature  in  Kelvin. 

The  sensible  heat  flux  is  relative  to  the  difference  between  air  potential  temperature  and  the 
surface  temperature  of  the  snow  or  snow-free  ice, 

F  =C  (0  -TKf). 

s  s  \  a  sf  ) 
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Cs  and  Ci  (below)  are  nonlinear  turbulent  heat  transfer  coefficients  described  in  Section  5.2. 1.1. 
Likewise,  the  latent  heat  flux  is  proportional  to  the  difference  between  Qa  and  the  surface 
saturation  specific  humidity  Qsf : 


F,  =  C,(Q,-Q„), 

Qsf  =  (yJpa)QM-chiTf), 

7  3 

where  qi  =1. 16378  x  10  kg/m  ,  q2  =  5897. 8K,  77  is  the  surface  temperature  in  Kelvin,  and  pa  is 
the  surface  air  density.  The  net  downward  heat  flux  from  the  ice  to  the  ocean  is  provided  by  [29]: 

Fbo,  =  ~Pwcwchu*  (r,  -Tf\  (6 !) 

where  pw  is  the  density  of  seawater,  cw  is  the  specific  heat  of  seawater,  c/,  =0.006  represents  a 
heat  transfer  coefficient,  n„  =  ■sj\fw\/pw  is  the  friction  velocity,  and  Tw  is  the  sea  surface 

temperature.  The  minimum  value  of  u„  depends  on  if  the  model  is  run  coupled;  lack  of  currents 
in  uncoupled  runs  means  there  was  insufficient  heat  available  to  melt  ice  in  the  standard 
formulation.  In  this  release  we  have  w*min  =  5xlCT3for  coupled  runs  and  5xlCT2  for  uncoupled 
runs. 


5.2.2. 5.2  New  Temperatures 

Given  the  temperatures  T"f‘ ,  T'"  and  T”  at  time  m,  we  solve  a  series  of  finite-difference  equations 

to  obtain  the  new  temperatures  at  time  m  +1.  Each  temperature  is  coupled  to  the  temperatures  of 
the  layers  directly  above  and  below  by  heat  conduction  terms  that  are  treated  implicitly.  For 
example,  the  rate  of  change  of  7T  depends  on  the  new  temperatures  in  layers  k- 1,  k,  and  k  + 1. 
Therefore,  we  have  a  set  of  equations  of  the  form 

Ax  =  b,  (62) 

where  A  is  a  tridiagonal  matrix,  x  is  a  column  vector  whose  components  are  the  unknown  new 
temperatures,  and  b  represents  another  column  vector.  Given  A  and  b,  x  can  be  computed  with  a 
standard  tridiagonal  solver. 

There  are  four  general  cases: 

(1)  Tsf<  0°C,  snow  present; 

(2)  Tsf=  0°C,  snow  present; 

(3)  Tsf<  0°C,  snow  absent;  and 
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(4)  Tsf= 0°C,  snow  absent. 

Case  1  has  one  equation  (the  top  row  of  the  matrix)  for  the  new  surface  temperature,  one 
equation  (the  second  row)  for  the  new  snow  temperature,  and  /V,  equations  (the  remaining  rows) 
for  the  new  ice  temperatures.  For  cases  2  and  4  the  equation  for  the  surface  temperature  is 
omitted,  which  is  held  at  0°C,  and  for  cases  3  and  4  the  snow  temperature  equation  is  omitted. 


The  rate  of  temperature  change  in  the  ice  interior  is  given  by  [31]: 


Pfi 


dT: 


dt  dz 


d  f.  8T^ 
ki — L 
dz 


—Uo  exp(-/qz)], 
dz 


(63) 


3 

where  pt  =  917  kg/m  represents  the  sea  ice  density  (assumed  to  be  unifonn),  Cj(T,S)  is  the 
specific  heat  of  sea  ice,  kj(T,S)  is  the  thermal  conductivity  of  sea  ice,  Iq  is  the  flux  of  penetrating 
solar  radiation,  attenuated  with  extinction  coefficient  k.  (see  previous  section),  and  z  represents 

the  vertical  coordinate,  set  to  be  positive  downward  with  z  =0  at  the  top  surface.  The  specific 
heat  of  sea  ice  is  given  to  an  excellent  approximation  by  [33], 


c,(r,S)  =  c„+^,  (64) 

5 

where  Co  =  2106  J/kg/deg  is  the  specific  heat  of  fresh  ice  at  0°C,  L0  =3.34  x  10  J/kg  is  the  latent 
heat  of  fusion  of  fresh  ice  at  0°C,  and  p  =0.054  deg/psu  represents  the  ratio  between  freezing 
temperature  and  the  salinity  of  brine.  Following  [44],  the  thermal  conductivity  is  shown  by 


k,  (T,S)  —  k0  + 


§S_ 

T  ’ 


(65) 


where  k0  =2.03  W/m/deg  is  the  conductivity  of  fresh  ice  and  />  =0. 13  W/m/psu  is  an  empirical 
constant.  The  analogous  equation  for  the  temperature  change  in  snow  is 


PsCs 


dT 


dt  dz 


d  ( ,  dT  ) 


(66) 


3 

where  ps  =  330  kg/m  is  the  snow  density  (also  assumed  uniform),  cs  =  c0  is  the  specific  heat  of 
snow,  and  ks  =  0.30  W/m/deg  denotes  the  thennal  conductivity  of  snow.  Penetrating  solar 
radiation  is  neglected  in  equation  (66)  because  the  majority  of  the  incoming  sunlight  is  absorbed 
near  the  top  surface  when  snow  exists. 


Now  Equations  (63)  and  (66)  are  converted  to  finite-difference  form.  The  resulting  equations  are 
second-order  accurate  in  space  (except  maybe  at  material  boundaries)  and  first-order  accurate  in 
time.  Before  writing  the  equations  in  full  the  finite-difference  expressions  are  provided  for  some 
of  the  terms. 
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First  consider  the  terms  on  the  left-hand  side  of  equations  (63)  and  (66).  The  time  derivatives  are 
written  as 


8T 


-T/W+l 


T 


m 


dt  At 


where  T  is  the  temperature  of  either  ice  or  snow.  The  specific  heat  of  ice  layer  k  is  approximated 
as 


:  cn  + 


LoMSa 

rpm  rpm-V 
1ik  1ik 


(67) 


which  guarantees  that  energy  is  conserved  during  a  change  in  temperature.  This  can  be  shown  by 
using  equation  (64)  to  integrate  c,  dT  from  T™  to  Tjk+1 ;  the  result  is  cik(T£+1  -71")  ,  where  c is 

given  by  equation  (67).  Unfortunately,  the  specific  heat  is  a  nonlinear  function  of  T"‘ 1 1 ,  the 
unknown  new  temperature.  A  set  of  linear  equations  is  retained,  however,  by  initially  guessing 
Tjk+l  =  T"k  and  then  iterating  the  solution,  updating  Tik+1  in  equation  (67)  with  each  iteration  until 
the  solution  converges. 


Next  to  consider  is  the  heat  diffusion  tenn,  the  first  term  on  the  right  side  of  equation  (63).  In  the 
ice  interior  (layers  2  to  Nj-  1)  this  term  is  discretized  as 


1  1 

1  / rpm+1  rpm+ 1\ 

Ki,k+ 1  ik  ~  1  i,k+ 1 ) 

V.  '  Sz  ) 

>  1 

A  ht 

> 

.s- 

1 _ 

(68) 


where  k %  represents  the  thermal  conductivity  at  the  upper  boundary  of  layer  k.  The 
approximation  in  equation  (68)  is  spatially  centered  and  second-order  accurate.  Similar 
expressions  can  be  written  for  heat  diffusion  in  the  top  and  bottom  ice  layers  and  the  snow  layer, 
as  shown  below.  Note  that  the  conduction  terms  are  handled  implicitly;  that  is,  they  depend  on 
the  temperatures  at  the  new  time  m+ 1.  The  resulting  scheme  is  first-order  accurate  in  time  and 
unconditionally  stable. 

Using  equation  (65),  k is  approximated  in  the  ice  interior  (at  the  upper  boundary  of  layers  2  to 
N,)  as 


,  *  -  -V,) 

kjk  —  kn  H - 

IK  U  rpm  ,  rpm 

1i,k-\  +  1ik 
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Because  the  conductivity  does  not  depend  as  much  on  temperature  as  does  the  specific  heat,  kik  is 
defined  in  terms  of  the  ice  temperatures  at  time  m.  Thus  the  conductivity  does  not  need  to  be 
updated  with  each  iteration.  At  the  bottom  surface  there  is 


:  K  + 


fisn 


The  conductivity  at  the  top  ice  surface,  ka ,  depends  on  whether  snow  exists.  If  there  is  no  snow, 
we  set 

kn=k0  +  —. 

ii  U  rnm 


(Defining  kn  in  terms  of  Tsf  is  avoided  since  then  it  would  be  undefined  for  Tsf  =0.)  If  snow  is 
present,  a  continuous  heat  flux  across  the  ice-snow  interface  is  assumed: 


k,  i 


rpm  rpm 
1i\  1int  _  K 

A/7,/2 


rpm  rpn 
1  int  1  s 

hJ2 


where  Tint  is  the  interface  temperature.  Solving  for  T™t ,  it  is  evident  that  this  heat  flux  is 
equivalent  to 

rpm  rpm 

k  ln  ~ls 
im(Ahi+hs)/2’ 

where  kint ,  the  equivalent  conductivity  at  the  interface,  is  defined  as 

k  _  WM  +K) 

M  Kkn+Akks  ’ 

Finally,  the  second  term  on  the  right  in  equation  (63)  is  taken  into  account.  From  equation  (60), 
the  fraction  of  the  penetrating  solar  radiation  Io  transmitted  through  layer  k  without  being 
absorbed  is 


rk  =  cxpi-k^kA/t ). 

Thus  the  flux  absorbed  in  layer  k  is 

Qk=I o(Vi-Tt)- 

The  flux  absorbed  per  unit  ice  thickness  is  Qi/Ahi,  the  desired  finite-difference  approximation  to 
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8 

-  —  [^o  exp(-K-,.z)]. 

dz 

Now  a  system  of  equations  for  the  new  temperatures  is  constructed.  (The  reader  uninterested  in 
algebraic  details  may  want  to  go  to  the  next  section.)  Beginning  at  the  surface  and  working  down 
for  case  1  (Tsf<  0°C  and  snow  present),  we  require 

F,=FC„  (69) 

where  Fct  is  the  conductive  flux  from  the  top  surface  to  the  ice  interior,  and  both  fluxes  are 
evaluated  at  time  m  +1.  Although  Fq  is  a  nonlinear  function  of  Tsfi  there  is  the  linear 
approximation 

Fm+l  =  p*  + 

1  0  1  0  T 

where  T*  is  the  surface  temperature  from  the  most  recent  iteration,  and  F0*  and  (dFJdTsf)  are 
functions  of  T*f  .  We  initialize  T*f  =  T"‘  and  update  it  with  each  iteration.  Thus  we  rewrite 
equation  (69)  as 


dF0 
dT , 


c Kr+l-T;f ) 
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Sf> 


F0  + 


(  A 
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\dTsf  J 


{TTx-rsf)=Ks{Tr-Tr) 


where  K  =  2kjhs .  Rearranging  terms,  we  get 


f  A 

dFn 
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(70) 


the  first  equation  in  the  set  of  equations  (62). 

Continuing  with  case  1,  we  write  the  equation  for  the  change  in  Ts: 


PsCs 


(Tm+I  -  Tn)  _  1 


At 


-[ks(t;+>-ts 


m+ 1 


)-Kw,(Tr-T;n\ 


where  Kint  =  2kinl/(Ah  +  h  ) .  In  tridiagonal  matrix  fonn  equation  (71)  becomes 


+  Ks)]Tsm+-TjsKintT” 


_  rjin 


(71) 


(72) 
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where  r/s  =  At/(pscshs) . 

The  ice  equations  for  the  top  layer,  the  interior  layers  (2  to  Nt  -  1),  and  the  bottom  layer, 
respectively,  are 


Pici 
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where  Kik  =  kikl Ah t  and  N  =  Nr  The  coefficients  yi  =3  and  }>2  =  -1/3  provide  one-sided  second- 

order  spatial  accuracy  at  the  bottom  surface;  they  are  derived  from  a  Taylor  series  expansion  of 
dT/dz  at  z  =  /?,.  Rearranging  terms,  we  have 


-VaKintTrl  +[l  +  7fl(^  +Kl2Wr =  Vi  +  7a&. 


m+\ 


im+l  _  rrim 


-VATZX  +  [!  +  %(*„ -pkK,Mir:Z  =  V+rjikQk, 


-^AKi,N-r2Ki,N+ +[1+7W(^  +r,^,  vuK 


/M+l 


?IiNKi.N+ 1  O')  +  /2  )7/  +  +  VinQn  ’ 


(73) 

(74) 

(75) 


where  rjik  =  At!{pfikAh/) . 

Next  consider  case  2  ( 7)/  =0°C  and  snow  present).  Since  Tsf  is  fixed,  there  is  no  surface  flux 
equation.  The  new  snow  temperature  is  given  by  equation  (71),  but  with  the  unknown 
TsJ+1  replaced  by  Tsf  =0°C.  In  matrix  form  we  have 

[1  +  17s(Kint  +  Ks)]Trl -dsK,Jn+l  =  VsKsTsf+V.  (76) 

The  ice  equations  for  case  2  are  the  same  as  for  case  1:  (73),  (74),  and  (75). 

For  case  3  ( Tsf<  0°C  and  snow  absent)  the  surface  temperature  equation  is  like  equation  (69),  but 
we  use  a  second-order  accurate  expression  for  dT/dz  at  z  =0: 
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^0  + 


rdF  Y 

-tt  0T1  - r;f)  =  kmt?'  - tt1) + r2(T£+l  - zr1)]. 

\dTsf) 


This  gives  the  matrix  equation 


f  A 

dF, 


\dTsf  J 


~Kn(ri+r2) 


Tsr'+nKnTrl+r2KnTr2+i  = 


f  A 
dFn 


ydTsfj 


T  -F 

1sf  1  0 


(77) 


The  equation  for  Tn  ty’in  case  3  is 


/ rjim+l  rr-im  \  -I 

pfn—  ~  — =Tr^[r1(T1-^+1)+^(T1-^r1)]-^2(^r1-^r1)+a- 

Ar  An. 


Rearranging  terms,  we  find 


-rt(7,+72)r  +[l  +  %(^2  +yAi)Ffl"+1  -%(^2 -r2Kn)Ttr  =  C  +  7„0,.  (78) 


Equation  (77)  includes  71™ +1  and  therefore  provides  an  unwanted  matrix  term  two  places  to  the 
right  of  the  main  diagonal.  This  term  is  eliminated  by  making  the  substitution 

7?j  — ^  c2R\  — CjT?, , 


where  Ri  is  the  first  matrix  row,  IF  is  the  second  row,  and  c]  =  y2 Kn  and  c2  =  -ijn(Kn  —y2Ka) 
are  the  coefficients  multiplying  77?+1  in  rows  1  and  2,  respectively.  The  other  ice  layer  equations 
for  case  3  are  equations  (74)  and  (75). 

Finally,  for  case  4  (Tsf=  0°C  and  snow  absent)  we  have  the  top  ice  layer  equation 


Pfi\ 


/ rjim+l  rrim\ 

vAl  ~1i\  ) 


i\  /  _ 


1 


At 


=  M  KMT*~T' 


)  +  n(T„ -T-”)]-Kn(T~" -T-")  +  Qv 


which  can  be  rewritten  as 


[  1  +  Vn  +  Kn  )]T™  +  rjn  (y2Kn -Kn)  T™  =  rjnKn(rl+y2)Tsf +T” +rjnQv 


The  remaining  ice  layer  equations  are  (74)  and  (75),  as  with  the  other  three  cases. 


(79) 
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This  completes  the  specification  of  the  matrix  equations  for  the  four  cases.  The  new  temperatures 
are  computed  using  a  tridiagonal  solver.  After  each  iteration  there  is  a  check  to  see  whether  the 
following  conditions  hold: 

1.  Tsf  <  0°C 

2.  The  change  in  Tsf  since  the  previous  iteration  is  less  than  a  prescribed  limit,  A Tmax. 

3.  F0>Fct.  (If  F0  <  Fct ,  ice  would  be  growing  at  the  top  surface,  which  is  not  allowed.) 

4.  The  rate  at  which  energy  is  added  to  the  ice  by  the  external  fluxes  equals  the  rate  at  which 
the  internal  ice  energy  is  changing,  to  within  a  prescribed  limit  A Fmax. 


The  convergence  rate  of  Tsf  is  also  checked.  If  Tsf  is  oscillating  and  failing  to  converge, 
temperatures  from  successive  iterations  are  averaged  to  improve  convergence.  When  all  these 
conditions  are  satisfied,  typically  within  two  to  four  iterations  for  AT  iax  ~  0.0 1C  and 

A Fmax  ~  0.01  W/m2 ,  the  computation  is  complete. 

5.2.2.53  Growth  and  Melting 

First  the  expressions  are  derived  for  the  enthalpy  q.  The  enthalpy  of  snow  (or  fresh  ice)  is  given 
by 

ds(T)  =  -ps(-cQT  +  L0 ). 


Sea  ice  enthalpy  is  more  complex,  due  to  brine  pockets  whose  salinity  varies  inversely  with 
temperature.  The  specific  heat  of  sea  ice,  shown  by  equation  (64),  includes  the  energy  needed  to 
warm  or  cool  ice  as  well  as  the  energy  used  to  freeze  or  melt  ice  adjacent  to  brine  pockets. 
Equation  (64)  can  be  integrated  to  give  the  energy  8e  required  to  raise  the  temperature  of  a  unit 
mass  of  sea  ice  of  salinity  S  from  T  to  T  : 


Se(T,T  )  =  c0(T  -T)  +  L0/jS 


n 

kT 


_n 

T's 


If  we  let  T  =Tm=  -pS ,  the  temperature  at  which  the  ice  is  completely  melted,  we  have 


Se(T,TJ  =  c0(Tm-T)  +  L0(  1 

V 


3 

Multiplying  by  pl  to  transfonn  the  units  from  J/kg  to  J/m  and  adding  a  term  for  the  energy 
required  to  raise  the  meltwater  temperature  to  0°C,  we  get  the  sea  ice  enthalpy: 
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q,(T,S)  =  -Pi 


(  j  \ 

Pi 

co(Tm  ~T)  +  L0 

1  _  m 

1  ml 

~cTm. 

l  T  J 

(80) 


Note  that  equation  (80)  is  a  quadratic  equation  in  T.  Given  the  layer  enthalpies,  the  temperatures 
can  be  computed  using  the  quadratic  formula: 

_  -b-  ^ lb 2  -4 ac 
2a  ’ 


where 


Cl  Cq  ? 


b  =  (cw-c0)Tm-  —  -L0, 
Pi 

C  =  L0Tm- 


The  other  root  is  unphysical. 

Melting  at  the  top  surface  is  demonstrated  by 


(F0-Fci)At  ijF0  >  Fct 
0  otherwise 

where  q  is  the  enthalpy  of  the  surface  ice  or  snow  layer  (recall  that  q 
thickness.  If  the  layer  melts  completely,  the  residual  flux  is  used  to 
energy  remaining  when  the  ice  and  snow  have  melted  is  added  to 
cannot  grow  at  the  top  surface,  but  new  snow  can  fall.  Snowfall 
thennodynamic  time  step. 

Growth  and  melting  at  the  bottom  ice  surface  are  given  by 

q8h  =  (Fcb  -Fbol)At,  (82) 

where  Fbot  is  given  by  equation  (61)  and  Fcb  is  the  conductive  heat  flux  at  the  bottom  surface: 

Fcb  =  ^[7ATiN-Tf)  +  y2(T,N_x-Tf)}. 

If  ice  is  melting  at  the  bottom  surface,  then  q  in  equation  (82)  is  the  enthalpy  of  the  bottom  ice 
layer.  If  ice  is  growing,  q  is  the  enthalpy  of  new  ice  with  temperature  7)  and  salinity  Smax.  This 
ice  is  added  to  the  bottom  layer. 


<  0)  and  Sh  is  the  change  in 
melt  the  layers  below.  Any 
the  ocean  mixed  layer.  Ice 
is  added  at  the  end  of  the 
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If  the  latent  heat  flux  is  negative  (i.e.,  latent  heat  is  transferred  from  the  ice  to  the  atmosphere), 
snow  or  snow-free  ice  sublimates  at  the  top  surface.  If  the  latent  heat  flux  is  positive, 
atmospheric  vapor  is  deposited  at  the  surface  as  snow  or  ice.  The  thickness  change  of  the  surface 
layer  is  given  by 


(pL-q)Sh  =  F,At, 


(83) 


where  p  is  the  density  of  the  surface  material  (snow  or  ice),  and  Lv  =  2.501x  106  J/kg  is  the  latent 

heat  of  vaporization  of  liquid  water  at  0°C.  Notice  that  pLv  is  close  to  an  order  of  magnitude 
larger  than  typical  values  of  q.  For  positive  latent  heat  fluxes,  the  deposited  snow  or  ice  is 
assumed  to  have  the  same  enthalpy  as  the  existing  surface  layer. 

After  growth  and  melting,  the  various  ice  layers  no  longer  have  the  same  thicknesses.  Therefore 
the  layer  interfaces  are  adjusted,  conserving  energy,  so  as  to  re-establish  layers  of  equal  thickness 
A  hi  =  hj/Ni.  This  is  accomplished  by  computing  the  overlap  qkm  of  each  new  layer  k  with  each  old 
layer  m: 

Vkm  =  min(zm  ,zk)~  ma x(zm_, ,  z,  , ), 


where  zm  and  zk  are  the  vertical  coordinates  of  the  old  and  new  layers,  respectively.  The 
enthalpies  of  the  new  layers  are 


1  «' 

CIk  .  ,  YjikmVm ' 

A  h 


At  the  end  of  the  time  step  there  is  a  check  of  whether  the  snow  is  deep  enough  to  be  partially 
below  freeboard  (i.e.,  below  the  surface  of  the  ocean).  Using  Archimedes’  principle,  the  base  of 
the  snow  is  at  freeboard  when 

M  +  PA  =  P,hr 


So  then  the  snow  base  is  below  freeboard  when 


h  =h- 


(Pu-PiVh 

Ps 


>0. 


In  this  case  the  snow  base  is  raised  to  freeboard  by  converting  some  snow  to  ice: 


Sh, 
8ht , 


~P,h* 

Pw 

PJL 

Pw 
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In  exceptional  cases  this  process  may  increase  the  ice  thickness  substantially.  Therefore  we 
postpone  snow-ice  conversions  until  after  the  remapping  in  thickness  space  (Section  5. 2. 2. 2), 
which  assumes  that  ice  growth  during  a  single  time  step  is  quite  small.  Lateral  melting  is 
accomplished  by  multiplying  the  state  variables  by  1  -rside,  where  rside  is  the  fraction  of  ice  melted 
laterally,  and  adjusting  the  ice  energy  and  fluxes  as  appropriate. 
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5.3  Primary  PIPS  3.0  FORTRAN  Routines 
5.3.1  PIPS  3.0  Modules 

Module  Description 

ice_albedo  Snow  and  ice  albedo  parameterization  and  aggregation. 

I/O:  None 

Calls:  icekindsmod,  icedomain 

Called  by:  icemodel,  ice  coupling,  inputdata,  mixed  layer,  absorbedsolar, 
scale_fluxes,  runtime_diags,  ice_write_hist 

I/O  Parameters:  None 

ice_atmo  Atmospheric  boundary  interface  (stability  based  flux  calculations). 

I/O:  None 

Calls:  ice  domain,  ice  constants,  ice  flux,  ice  state 

Called  by:  None 

I/O  Parameters:  None 

ice_calendar  Calendar  routines  for  managing  time. 

I/O:  None 

Calls:  iceconstants 

Called  by:  evp_prep,  dumpfile,  icecdf,  icemodel,  ice  coupling,  ice  diagnostics, 
ice  flux  in,  iceitdlinear,  ice  mechred,  ice  thenn  vertical,  icethennitd, 
ice_transport_remap,  ice_write_hist,  init  hist,  input_data,  init_evp,  mixed_layer, 
mpdata,  restartfile,  zapsmallareas 

I/O  Parameters:  None 

ice_constants  This  module  defines  a  variety  of  physical  and  numerical  constants  used 

throughout  the  ice  model. 

I/O:  None 

Calls:  If  coupled,  shrconstmod.  Otherwise,  ice  kinds  mod,  ice  domain 
Called  by:  albedos,  global_gather,  global_scatter,  ice_atmo,  icecdf,  ice_calendar, 
ice_coupling,  ice_diagnostics,  ice_dyn_evp,  ice_flux,  ice_flux_in,  ice_grid, 
ice  itd,  ice  itd  linear,  ice  mechred,  ice  ocean,  ice  scaling,  ice  thenn  itd, 
ice  thenn  vertical,  ice  timers,  icetransportmpdata,  icetransportremap, 
ice_write_hist,  init_flux,  init  hist,  init_state 
I/O  Parameters:  None 

ice_coupling  Message  passing  to  and  from  the  coupler. 

I/O:  None 

Calls:  ice  kinds  mod,  icemodelsize,  ice  constants,  ice  grid,  ice  state, 

ice_flux,  ice_albedo,  ice_mpi_intemal,  ice_timers,  ice_fileunits,  ice_work  (only: 

worka,  workjl).  If  coupled,  shr_sys_mod  (only:  shr_sys  Jlush), 

cpl  contract  mod,  cpl  interface  mod,  cpl  fields  mod 

Called  by:  icemodel,  dumpfile,  input  data,  restartfile,  setup  mpi 

I/O  Parameters:  None 

ice_diagnostics  Diagnostic  infonnation  output  during  a  model  run. 

I/O:  None 

Calls:  ice_domain,  ice_constants,  ice_calendar,  ice_fileunits,  ice_work  (only: 
work_g2,  work  ll,  work  12,  worka,  workb) 
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Module 


ice  domain 


ice_dyn_evp 


ice  exit 


ice Jileunits 


ice Jlux 


Description 

Called  by:  icemodel,  debug  ice,  icethennitd,  icethennvertical,  inputdata, 
readclimdata,  read  data 

I/O  Parameters:  None 

Sets  array  sizes  for  the  local  subdomain  and  related  parallel  processing 
infonnation.  This  code  was  originally  based  on  domain. F  in  the  POP  model. 

I/O:  None 

Calls:  icekindsmod,  icemodelsize 

Called  by:  abort_ice,  aggregate,  icemodel,  ice_atmo,  ice_dyn_evp,  ice_flux_in, 
ice_grid,  ice_history,  ice_init,  ice_itd_linear,  ice_mechred,  ice_read_write, 
ice_scaling,  ice_therm_itd,  ice_thenn_vertical,  ice_timer_print, 
icetransportmpdata,  icetransportremap 

I/O  Parameters:  None 

This  is  the  elastic -viscous-plastic  sea  ice  dynamics  model.  It  computes  ice 
velocity  and  defonnation. 

References: 

[18],  [17],  [19]  and  [20]. 

I/O:  None 

Calls:  ice  kinds  mod,  icedomain,  ice  grid,  ice  constants,  ice  state,  ice  work 
(only:  worka,  workb) 

Called  by:  icemodel,  dumpfile,  ice_write_hist,  input_data,  restartfde 

I/O  Parameters:  None 

Exits  the  model.  Logically,  this  routine  should  be  used  for  "normal"  exit  of  the 
ice  model,  but  there  would  be  only  one  such  call,  and  creating  a  subroutine  here 
to  accomplish  that  causes  circular  dependencies  due  to  the  coupler  exit  strategy. 
I/O:  None 

Calls:  icekindsmod 
Called  by:  None 
I/O  Parameters:  None 

Defines  unit  numbers  for  files  opened  for  reading  or  writing. 

I/O:  None 

Calls:  icekindsmod 

Called  by:  abort  ice,  calendar,  icemodel,  icecoupling,  ice  diagnostics, 
ice_flux_in,  ice_global_real_mimnax,  ice_grid,  ice_history,  ice_itd, 
ice_itd_linear,  ice_mechred,  ice_read_write,  ice_thenn_vertical,  ice_timer_print, 
ice  transport  mpdata,  ice  transport  remap,  init  evp 

I/O  Parameters:  None 

Flux  variable  declarations.  These  include  fields  sent  from  the  coupler  ("in"),  sent 
to  the  coupler  ("out"),  written  to  diagnostic  history  files  ("diagnostic"),  and  used 
internally  ("internal"). 

I/O:  None 

Calls:  ice  kinds  mod,  ice  domain,  ice  constants 

Called  by:  aggregate,  check_state,  dumpfile,  evp_finish,  evp_prep,  ice_atmo, 
ice  coupling,  ice  flux  in,  ice  scaling,  icethermitd,  ice  thenn  vertical, 
ice_write_hist,  init_evp,  init_flux,  init  hist,  init_state,  mixed  layer,  print_state, 
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Module 


ice Jluxjn 


ice_grid 


ice_histojy 


ice  init 


ice  itd 


ice  itd  linear 


Description 

ridge  ice,  runtimediags,  shift  ice,  stepu,  transportmpdata,  zapsmallareas 

I/O  Parameters:  None 

Reads  and  interpolates  forcing  data  for  atmospheric  and  oceanic  quantities. 

I/O:  None 

Calls:  icekindsmod,  icedomain,  ice  constants,  iceflux,  icecalendar 
ice_read_write,  ice_fdeunits 
Called  by:  icemodel,  inputdata 

I/O  Parameters:  None 

This  module  contains  spatial  grids,  masks,  and  boundary  conditions. 

I/O:  None 

Calls:  ice  kinds  mod,  ice  constants,  ice  domain,  ice  fl leunits, 
ice_mpi_intemal,  ice_work  (only:  work_gl,  work_g2,  work  ll,  worka) 

Called  by:  aggregate,  atmo_boundary_layer,  bound_aggregate,  bound_state, 
icemodel,  debug  ice,  dumpfile,  icecdf,  icemechred,  ice  scaling,  icethennitd, 
ice  therm  vertical,  icetransportmpdata,  icetransportremap,  ice  write  hist, 
init_state,  mixed_layer,  rebin,  reduce_area,  re  start  file 
I/O  Parameters:  None 

Reads/Writes  ice  model  history  and  restart  files.  Output  files  include  netCDF 
data  and  Fortran  unformatted  dumps. 

I/O:  None 

Calls:  ice  kinds  mod,  ice  domain,  ice  read  write,  ice  fileunits,  ice  work 
(only:  work_gl,  work_gr,  worka) 

Called  by:  icemodel,  input  data 
I/O  Parameters:  None 
Parameter  and  variable  initializations. 

I/O:  None 
Calls:  icedomain 

Called  by:  icemodel,  ice  transport  mpdata 

I/O  Parameters:  None 

Initializes  and  redistributes  ice  in  the  ITD.  Ice_itd  contains  routines  to  initialize 
the  ice  thickness  distribution  and  utilities  to  redistribute  ice  among  categories. 
These  routines  are  not  specific  to  a  particular  numerical  implementation. 
References:  [7],  [6]. 

I/O:  None 

Calls:  abort  ice,  ice  kinds  mod,  icemodelsize,  ice  constants,  ice  state, 
ice_fileunits, 

Called  by:  icemodel,  debug  ice,  ice  itd  bnear,  ice  mechred,  ice  thenn  itd, 

ice  therm  vertical,  ice  transport  remap,  init  thermo  vertical,  thenno  itd, 

transportmpdata 

I/O  Parameters:  None 

Linear  remapping  scheme  for  the  ITD. 

I/O:  None 

Calls:  ice  model  size  ice  kinds  mod,  ice  domain,  ice  constants  ice  state, 
ice  itd,  ice  calendar,  ice  fileunits 
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Module 

icejdndsjnod 

icejnechred 

icemodel 

ice_model_size 

ice  mpi internal 

ice  ocean 


Description 

Called  by:  icemodel,  thennoitd 

I/O  Parameters:  None 

Defines  variable  precision  for  all  common  data  types. 

I/O:  None 
Calls:  None 

Called  by:  icemodel,  debug_ice,  ice_albedo,  ice_constants,  ice_coupling, 
ice_domain,  ice_exit,  ice_dyn_evp,  ice_fileunits,  ice_flux,  ice_flux_in,  ice_grid, 
ice  history,  iceitd,  iceitdlinear,  icempiintemal,  ice  ocean,  ice  scaling, 
ice_state,  ice_thenn_itd,  ice_therm_vertical,  ice_timers,  ice_transport_remap, 
ice_work,  print_state 
I/O  Parameters:  None 

Ice  mechanical  redistribution  (ridging)  and  strength  computations. 

I/O:  None 

Calls:  ice_model_size,  ice_constants,  ice_state,  ice_itd,  ice_grid,  ice_fdeunits, 
ice_domain,  ice_calendar,  ice_work  (only:  worka,  workb ) 

Called  by:  icemodel 
I/O  Parameters:  None 

Main  driver  routine  for  PIPS  3.0.  Initializes  and  steps  through  the  model. 

I/O:  None 

Calls:  ice  albedo,  ice  calendar,  icecoupling,  ice  diagnostics,  icedomain, 

ice_dyn_evp,  ice_fdeunits,  ice_flux_in,  ice_grid,  ice_history,  ice_init,  ice_itd, 

ice  itd  linear,  icekindsmod,  icemechred,  ice  mpi  intemal,  ice  ocean, 

ice_scaling,  ice_therm_vertical,  ice_thenn_itd,  ice_timers,  ice_transport_mpdata, 

icetransportremap 

Called  by:  None 

I/O  Parameters:  None 

Defines  the  global  domain  size  and  number  of  categories  and  layers. 

The  code  is  based  on  model_size.F  in  the  POP  model. 

I/O:  None 

Calls:  icekindsmod 

Called  by:  columngrid,  dumpfile,  global  gather,  global  scatter,  icecdf, 
ice  coupling,  ice  domain,  ice  itd,  ice  itd  linear,  ice  mechred  ice  read  write, 
ice  state,  icethennitd,  ice  thenn  vertical,  icetransportmpdata, 
ice_transport_remap,  init_state,  print_state,  rectgrid,  restartfile,  runtime_diags, 
tlatlon 

I/O  Parameters:  None 

Parameters  and  common  blocks  for  MPI  parallelization  internal  to  the  ice  model. 
I/O:  None 

Calls:  ice  kinds  mod,  ice  domain 

Called  by:  abort_ice,  conserved_sums,  debug_ice,  icecdf,  icemodel, 
ice_coupling,  ice_grid,  ice_read_write,  ice_timer_print,  init_diags, 
init_mass_diags,  restartfile,  runtime_diags,  setup_mpi 

I/O  Parameters:  None 

Ocean  mixed  layer  calculation  (internal  to  sea  ice  model).  It  allows  heat  storage 
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Module  Description 

in  the  ocean  for  uncoupled  runs. 

I/O:  None 

Calls:  icekindsmod,  iceconstants 
Called  by:  icemodel,  thermovertical 

I/O  Parameters:  None 

Routines  for  opening,  reading  and  writing  external  files. 

I/O:  None 

Calls:  icemodelsize,  icedomain,  icempiintemal,  ice  fileunits,  ice  work 
(only:  work_gl,  work_gr ) 

Called  by:  ice_flux_in,  ice_history,  pipsgrid,  popgrid 
I/O  Parameters:  None 
Scales  ice  fluxes  by  ice  area. 

I/O:  None 

Calls:  ice  domain,  ice  kinds  mod,  ice  constants,  ice  state,  ice  flux,  ice  grid 
(only:  tmask ) 

Called  by:  icemodel 
I/O  Parameters:  None 

Primary  state  variables  in  various  configurations. 


The  primary  state  variable  names  are: 


For  Each  Category 

Aggregated  Over  Units 

Categories 

aicen(i,j,n) 

aice(i,j) 

— 

vicen(i,j,n) 

vice(ij) 

m 

vsnon(i,j,n) 

vsno(i,j) 

m 

eicen(i,j,k) 

eice(i,j) 

J/m2 

esnon(i,j,n) 

esno(i,j) 

J/m2 

Tsfcn(i,j,n) 

Tsfc(i,i) 

deg 

Area  is  dimensionless  because  aice  is  the  fractional  area  (normalized  so  that  the 
sum  over  all  categories,  including  open  water,  is  1.0).  That  explains  why 
vice/vsno  has  units  of  m  instead  of  m3,  and  eice/esno  have  units  of  J/m 2  instead  of 
J. 

I/O:  None 

Calls:  ice  kinds  mod,  ice  model  size,  ice  domain 

Called  by:  albedos,  dumpfile,  ice_atmo,  ice_coupling,  ice_dyn_evp,  ice_itd, 

ice_itd_linear,  ice_mechred,  ice_scaling,  ice_therm_vertical,  ice_therm_itd, 

icetransportremap,  ice  write  hist,  init  diagnostics,  init  flux  atm, 

init_mass_diags,  init_state,  merge_fluxes,  mixed_layer,  mpdata,  prepare_forcing, 

print  state,  restartfile,  runtime  diags,  transportmpdata 

I/O  Parameters:  None 

icethermitd  Thermodynamic  calculations  after  call  to  coupler,  mostly  related  to  ITD:  ice 

thickness  redistribution,  lateral  growth  and  melting,  and  freeboard  adjustment. 


ice  read  write 


ice_scaling 


ice  state 
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NOTE:  The  thermodynamic  calculation  is  split  in  two  for  load  balancing. 

First  icejhermvertical  computes  vertical  growth  rates  and  coupler  fluxes.  Then 
icejhermjtd  does  thennodynamic  calculations  not  needed  for  coupling. 

I/O:  None 

Calls:  icekindsmod,  icemodelsize,  ice  constants,  icedomain,  ice  state, 
ice_flux,  ice_diagnostics,  ice_calendar,  ice_grid,  ice_itd 
Called  by:  icemodel 

I/O  Parameters:  real  hicen-  ice  thickness  (m) 

icejhermjertical  Thermodynamic  calculations  before  calls  to  the  coupler.  This  module  updates  the 

ice  and  snow  internal  temperatures  and  computes  thennodynamic  growth  rates 
and  atmospheric  fluxes. 

References:  [7] 

NOTE:  The  thermodynamic  calculation  is  split  in  two  for  load  balancing. 

First  icejhermjertical  computes  vertical  growth  rates  and  coupler  fluxes.  Then 
icejhermjtd  does  thennodynamic  calculations  not  needed  for  coupling. 

I/O:  None 

Calls:  ice  model  size,  ice  kinds  mod,  ice  domain,  ice  fdeunits, 
ice_constants,  ice_calendar,  ice_grid,  ice_state,  ice_flux,  ice_itd,  ice_diagnostics 
(only:  print state) 

Called  by:  icemodel,  thennoitd 

I/O  Parameters:  None 

ice  timers  Timing  routines  for  the  PIPS  3.0  model. 

I/O:  None 

Calls:  ice  kinds  mod,  ice  constants 

Called  by:  bound  ijn,  evp,  icemodel,  decoupling,  icetransportremap, 
ridge  ice,  thenno  itd,  thermo  vertical,  transportmpdata 

I/O  Parameters:  None 

icejransportjnpdata  Calculates  horizontal  advection  using  MPDATA  (Multidimensional 

Positive  Definite  Advection  Transport  Algorithm). 

I/O:  None 

Calls:  ice  model  size,  ice  domain,  ice  constants,  ice  grid,  ice  fdeunits,  ice  init 
(only:  advection ),  ice_work  (only:  work 11,  work J 2) 

Called  by:  icemodel 

I/O  Parameters:  None 

ice  transport  remap  Transports  quantities  using  the  second-order  conservative  remapping  scheme 

developed  by  John  Dukowicz  and  John  Baumgardner  (DB)  and  modified  for  sea 
ice  by  William  Lipscomb  and  Elizabeth  Hunke. 

I/O:  None 

Calls:  ice  calendar,  ice  constants,  ice  domain,  ice  fdeunits,  ice  grid,  ice  itd, 
ice_kinds_mod,  icemodelsize,  ice_state,  ice_timers,  ice_work  (only:  work  11) 
Called  by:  icemodel 
I/O  Parameters:  None 

icejvork  The  intent  of  this  subroutine  is  to  have  one  global  work  array  available  all  the 

time,  and  another  available  that  can  be  allocated  when  necessary.  Globally 
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accessible,  local  (i.e.,  on-processor)  work  arrays  are  also  available  to  conserve 
memory.  These  arrays  should  be  used  only  within  a  single  subroutine. 

I/O:  None 

Calls:  ice  kinds  mod,  icedomain 

Called  by:  ice_coupling,  ice_diagnostics,  ice_dyn_evp,  ice_grid,  ice_history, 
icemechred,  ice  read  write,  icetransportmpdata,  icetransportremap, 
sss_clim,  thenno_vertical,  shift_ice 

I/O  Parameters:  None 


5. 3. 2  PIPS  3. 0  Subroutines 


Subroutine  Description 

abort  ice  This  routine  aborts  the  ice  model  and  prints  an  error 

message. 

I/O:  stdout 

Calls:  ice  domain,  icc  fl leun its,  icempiintemal.  If 
coupled,  shr  sys  mod 

Called  by:  aggregate_area,  check_monotonicity, 

columngrid,  columnconservationcheck, 

conservation  check  vthenno,  global  conservation, 

init  grid,  init  hist,  init  itd,  initremap,  init_vertical_profile, 

input  data,  mpdata,  ridge  ice,  ridge  shift,  setup  mpi 

shift_ice,  temperature_changes,  update_fields, 

zapsmallareas 

I/O  Parameters:  character  error_message  (input) 

absorbed  solar  Computes  solar  radiation  absorbed  in  ice  and  penetrating  to 

the  ocean. 

I/O:  None 

Calls:  icealbedo 

Called  by:  temperature_changes 

I/O  Parameters:  integer  n-  thickness  category  index 

integer  icells-  no.  of  cells  with  aicen  >  puny 

integer  indxi,  indxj-  compressed  indices  for  cells  with  aicen 

>  puny 

real  hlyr-  ice  layer  thickness 
real  hsn-  snow  thickness  (m) 

real  fswsfc-  SW  absorbed  at  ice/snow  surface  (W/m“) 
real  fswint-  SW  absorbed  in  ice  interior,  below  surface 
(W/m2) 
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add  new  ice 


add  new  snow 


aggregate 


aggregate _area 


Description 

real  fswthrun-  SW  through  ice  to  ocean  (W/m  ) 
real  labs-  SW  absorbed  in  particular  layer  (W/m2) 

Given  the  volume  of  new  ice  grown  in  open  water,  this 
subroutine  computes  its  area  and  thickness  and  adds  it  to  the 
appropriate  category  or  categories. 

NOTE:  Usually  all  the  new  ice  is  added  to  category  1.  An 
exception  is  made  if  there  is  no  open  water  or  if  the  new  ice 
is  too  thick  for  category  1,  in  which  case  ice  is  distributed 
evenly  over  the  entire  cell.  Subroutine  rebin  should  be 
called  in  case  the  ice  thickness  lies  outside  category  bounds 
after  new  ice  fonnation. 

I/O:  None 

Calls:  columnconservationcheck,  columnsum 
Called  by:  thennoitd 
I/O  Parameters:  None 

Adds  new  snow  at  top  surface. 

I/O:  None 
Calls:  None 

Called  by:  thennovertical 

I/O  Parameters:  integer  n  -  thickness  category  index 
integer  icells-  number  of  cells  with  aicen  >  puny 
integer  indxi,  indxj-  compressed  indices  for  cells  with  aicen 
>  puny 

real  Tsf-  ice/snow  surface  temperature,  Tsfcn 

real  hsn-  snow  thickness  (m) 

real  qsn-  snow  enthalpy 

real  hsn_new-  thickness  of  new  snow  (m) 

Aggregates  ice  state  variables  over  thickness  categories. 

I/O:  None 

Calls:  ice_domain,  ice_flux  (only:  Tf),  ice_grid 
Called  by:  icemodel,  init  state,  restartfile,  thenno  itd 

I/O  Parameters:  None 

Aggregates  ice  area  (but  not  other  state  variables)  over 
thickness  categories. 

I/O:  stdout 
Calls:  abort  ice 
Called  by:  linearitd 
I/O  Parameters:  None 
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albedos 


asum  ridging 


atmo _boundary Jay  er 


Description 

Computes  snow/ice  albedos  and  aggregate.  Ice  albedo  is  zero 
if  no  ice  is  present. 

Calls:  ice_constants,  ice_grid,  ice_state 
Called  by:  icemodel 
I/O  Parameters:  None 

Finds  the  total  area  of  ice  plus  open  water  in  each  grid  cell. 
This  is  similar  to  the  aggregate _area  subroutine  except  that 
the  total  area  can  be  greater  than  1,  so  the  open  water  area  is 
included  in  the  sum  instead  of  being  computed  as  a  residual. 
I/O:  None 
Calls:  None 

Called  by:  ridge_ice,  ridge_prep 

I/O  Parameters:  None 

Computes  coefficients  for  atmosphere-ice  fluxes,  stress  and 
2  meter  reference  temperature. 

NOTE: 

(1)  all  fluxes  are  positive  downward. 

(2)  net  heat  flux  =fswabs  +  jlwup  +  ( 1  -emissivity)/7>v<fn  + 
fsh  +  flh. 

(3)  here,  tstar  =  (WT)/U*,  and  qstar  =  (WQ)/U*. 

(4)  wind  speeds  should  all  be  above  a  minimum  speed  (eg. 
1.0  m/s). 

I/O:  None 
Calls:  ice_grid 

Called  by:  mixed  layer,  thennovertical 
Assumptions:  The  saturation  humidity  of  air  at  T(K): 
qsat(T)  (kg/m3). 

I/O  Parameters:  integer  n-  thickness  category  index 

character  sfctype-  ice  or  ocean 

real  Tsf-  surface  temp  of  ice  or  ocean 

real  strx-  x  surface  stress 

real  stry-  y  surface  stress 

real  Trf-  reference  height  temp  (K) 

real  Qrf-  reference  height  specific  humidity  (kg/kg) 

real  delt-  potential  T  difference  (K) 

real  delq-  humidity  difference  (kg/kg) 
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bound 


bound aggregate 


boundjjn 


bound  nan- 


bound  narr  ne 


Description 

Fills  ghost  cells  with  boundary  infonnation. 

I/O:  None 
Calls:  boundjjn 

Called  by:  bound_aggregate,  evp_prep,  from_coupler, 
init  evp,  initgrid,  initremap,  init  state,  local  max  min, 
makemask,  mpdata,  restartfde,  transport_remap,  t2ugrid, 
tlatlon,  tocoupler,  u2tgrid 
I/O  Parameters:  real  workl 

Bound  calls  for  aggregate  ice  state.  Gets  ghost  cell  values 
for  aggregate  ice  state  variables. 

NOTE:  This  subroutine  is  called  only  at  initialization. 

I/O:  None 

Calls:  ice  grid,  bound 
Called  by:  init_state,  restartfde 

I/O  Parameters:  None 

Periodic/Neumann  conditions  for  global  domain  boundaries. 
Assumptions:  A  “single”  row  of  ghost  cells  ( num-ghost - 
cells=l );  workl  array  has  fonn  (i-indexj -index,  number- 
arrays). 

I/O:  None 
Calls:  ice  timers 

Called  by:  bound,  bound  narr,  bound  narr  ne,  bound  sw 
I/O  Parameters:  integer  nd,  real  workl,  logical  north,  south, 
east,  west 

Fills  neighboring  ghost  cells  with  boundary  infonnation  and 
fills  several  arrays  at  once  (for  perfonnance). 

I/O:  None 
Calls:  boundjjn 

Called  by:  bound  state,  mpdata,  transportremap, 
transportmpdata 

I/O  Parameters:  integer  narrays,  real  workl 

Fills  north  and  east  ghost  cells  with  boundary  infonnation 
and  fills  several  anays  at  once  (for  perfonnance). 

I/O:  None 
Calls:  boundjjn 
Called  by:  stepu 

I/O  Parameters:  integer  nanays,  real  workl 
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bound  state 


bound  sw 


calendar 


checkmonotonicity 


check  state 


Description 

Bound  calls  for  ice  state  variables.  Gets  ghost  cell  values  for 
ice  state  variables  in  each  thickness  category. 

I/O:  None 

Calls:  ice  grid,  bound  narr 
Called  by:  restartfile,  transport_remap 

I/O  Parameters:  None 

Fills  south  and  west  ghost  cells  with  boundary  information. 

I/O:  None 

Calls:  bound  ijn 

Called  by:  evp,  transport_remap 

I/O  Parameters:  real  workl 

Determines  the  date  at  the  end  of  the  time  step. 

I/O:  stdout 

Calls:  ice _ H  leunits 

Called  by:  icemodel,  init  cpl 

I/O  Parameters:  real  ttime-  time  variable 

At  each  grid  point,  this  subroutine  assures  that  the  new  tracer 
values  fall  between  the  local  maximum  and  minimum  values 
before  transport. 

I/O:  stdout 

Calls:  abort  ice 

Called  by:  transportremap 

I/O  Parameters:  real  aim-  new  ice  area 

real  tnn-  new  tracers 

real  tmin-  local  min  tracer 

real  tmax-  local  max  tracer 

This  subroutine  requires  certain  fields  to  be  monotone. 
NOTE:  This  should  not  be  necessary  if  all  is  well,  but  it  is 
best  to  keep  going.  The  model  will  not  conserve  energy  and 
water  if  fields  are  zeroed  here. 

I/O:  None 

Calls:  iceflux 

Called  by:  transportmpdata 

I/O  Parameters:  None 
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column  conservation  check 


column  sum 


columngrid 

complete _getflux_ocn 


Description 

For  each  physical  grid  cell,  the  routine  checks  that  initial  and 
final  values  of  a  conserved  field  are  equal  to  within  a  small 
value. 

I/O:  stdout 
Calls:  abortice 

Called  by:  add_new_ice,  linear_itd,  ridge_shift 

I/O  Parameters:  real  xl-  initial  field 

real  x2-  final  field 

real  maxerr-  max  allowed  error 

character  fielded-  field  identifier 

For  each  grid  cell,  the  subroutine  sums  the  field  over  all  ice 
categories. 

I/O:  None 
Calls:  None 

Called  by:  add_new_ice,  linear_itd,  ridge_shift 
I/O  Parameters:  integer  nsum-  no.  of  categories/layers 
real  xin-  input  field 
real  xout-  output  field 

Constructs  column  grid  and  mask. 

I/O:  stdout 

Calls:  global  scatter,  icemodelsize,  abort  ice 
Called  by:  init_grid 
I/O  Parameters:  None 

Computes  remaining  ocean  forcing  fields. 

I/O:  None 
Calls:  None 

Called  by:  sss_clim,  sss_sst_restore 

I/O  Parameters:  None 
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conductivity 


conservation  check  vthermo 


Description 

Computes  thermal  conductivity  at  interfaces  (held  fixed 
during  the  subsequent  iteration). 

NOTE:  Ice  conductivity  must  be  >=  kimin 
I/O:  None 
Calls:  None 

Called  by:  temperature_changes 

I/O  Parameters:  integer  icells-  no.  of  cells  with  aicen  > 

puny 

integer  indxi,  indxj-  compressed  indices  for  cells  with  aicen 

>  puny 

real  hlyr-  ice  layer  thickness 

real  hsn-  snow  layer  thickness 

real  Tbot-  ice  bottom  surface  temp  (°C) 

real  Tin-  internal  ice  layer  temperatures 

real  khi-  ki/hlyr 

real  khs-  ksno/hsn 

Checks  for  energy  conservation  by  comparing  the  change  in 
energy  to  the  net  energy  input. 

I/O:  stdout 

Calls:  print_state,  abort_ice 
Called  by:  thennovertical 

I/O  Parameters:  integer  n-  thickness  category  index 
(diagnostic  only) 

integer  icells-  number  of  cells  with  aicen  >  puny 
integer  indxi,  indxj-  compressed  indices  for  cells  with  aicen 

>  puny 

real  fsurf  -  net  flux  to  top  surface,  not  including  fcondtop 
real  flatn  -  surface  downward  latent  heat  (W/m2) 
real  fhnetn-  fbot,  corrected  for  any  surplus  energy 
real  fswint-  SW  absorbed  in  ice  interior,  below  surface 
(W/m2) 

real  einit  -  initial  energy  of  melting  (J/m  ) 
real  efinal-  final  energy  of  melting  (J/m") 
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conserved  sums 


construct Jields 


debug Jce 


departure _points 


Description 

Computes  global  sums  of  conserved  variables  over  the 
physical  grid. 

I/O:  None 

Calls:  icempiinternal,  iceglobalrealsum 

Called  by:  transport_remap 

I/O  Parameters:  real  aim-  mean  ice  area 

real  asum-  global  sum  of  area 

real  tnn-  mean  tracer 

real  atsum-  global  sum  of  area*tracer 

Constructs  fields  of  ice  area  and  tracers. 

I/O:  None 

Calls:  limitedgradient 

Called  by:  transport_remap 

I/O  Parameters:  real  aim-  mean  ice  area 

real  aimask-  =  1.  if  ice  is  present,  =  0.  otherwise 

real  tnn-  mean  tracer 

real  trmask-  =  1 .  if  tracer  is  present,  =0.  otherwise, 
real  aic-  ice  area  at  geometric  center  of  cell 
real  aix,  aiy-  limited  derivative  of  ice  area  with  respect  to  x 
and  y 

real  trc-  tracer  at  geometric  center  of  cell 

real  trx,  try-  limited  derivative  of  tracer  with  respect  to  x  and 

y 

Wrapper  for  the  printstate  debugging  routine.  It  is  useful 
for  debugging  in  the  main  driver. 

I/O:  None 

Calls:  icekindsmod,  ice  itd,  ice  diagnostics, 
icempiinternal,  ice  grid,  print  state 
Called  by:  None 

I/O  Parameters:  character  (char  len),  intent(in) ::  plabeld 

Given  velocity  fields  on  cell  corners,  this  subroutine 
computes  departure  points  of  trajectories  using  a  midpoint 
approximation. 

I/O:  stdout 
Calls:  None 

Called  by:  transport_remap 

I/O  Parameters:  real  dpx,  dpy-  x  andy  coordinates  of 
departure  points  at  cell  comers 
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dumpfile 


end  run 


evp 


evp Jinish 


evp _prep 


exit_coupler 


Description 

Dumps  all  values  needed  for  a  restart. 

I/O:  stdout,  write  filename,  nudump 
Calls:  icemodelsize,  ice  flux,  ice  grid,  ice  calendar 
ice_state,  ice_dyn_evp,  ice_ocean  (only:  oceanmixedjce), 
ice_open,  ice_coupling  (only:  l_coup/ed),  ice_write 
Called  by:  icemodel 
I/O  Parameters:  None 

Ends  run  by  calling  mpi  finalize. 

I/O:  None 

Calls:  mpi_finalize 

Called  by:  icemodel,  abort  ice 

I/O  Parameters:  None 

Elastic-viscous-plastic  dynamics  driver. 

I/O:  None 

Calls:  ice-timers,  ice_timer_start,  evp_prep,  stepu,  stress, 

bound  sw,  evp  finish,  icetimerstop 

Called  by:  icemodel 

I/O  Parameters:  integer  kstmgth-  input 

Calculates  ice-ocean  stress. 

I/O:  None 

Calls:  ice_flux,  iceumask,  u2tgrid 

Called  by:  evp 

I/O  Parameters:  None 

Computes  quantities  needed  in  the  stress  tensor  (sigma)  and 
momentum  (u)  equations,  but  the  following  quantities  do  not 
change  during  the  thermodynamics/transport  time  step: 
wind  stress  shift  to  U  grid,  ice  mass  and  ice  extent  masks, 
pressure  (strength),  and  part  of  the  forcing  stresses,  initializes 
ice  velocity  for  new  points  to  ocean  surface  current. 

I/O:  None 

Calls:  bound,  tmask,  t2ugrid,  tougrid,  iceumask, 
ice_strength,  ice_flux,  ice_calendar,  ice_mechred 
Called  by:  evp 

I/O  Parameters:  integer  kstmgth-  input 

Exits  from  coupled/MPI  environment. 

I/O:  stdout 

Calls:  cpl  interface  fuialize,  mpiabort 
Called  by:  icemodel 

I/O  Parameters:  None 
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file_year 


fitjine 


fluxjntegrals 


Description 

Constructs  the  correct  name  of  the  atmospheric  data  file  to 
be  read,  given  the  year  and  assuming  the  naming  convention 
of  filenames  ending  with  'yyyy.daf . 

I/O:  stdout  (a,  14.4,  a) 

Calls:  None 

Called  by:  ncar  files,  read  data 
I/O  Parameters:  character  data  file 

Fits  g(h)  with  a  line,  satisfying  area  and  volume  constraints. 
To  reduce  round  off  errors  caused  by  large  values  of  gO  and 
gl,  it  computes  g(eta),  where  eta  =  h  -  hL,  and  hL  is  the  left 
boundary. 

I/O:  None 
Calls:  None 
Called  by:  linearitd 

I/O  Parameters:  integer  n-  category  index 
real  HbL,  HbR-  left  and  right  cat  boundaries 
real  hice-  ice  thickness 

real  gO,  gl-  coefficients  in  linear  equation  for  g(eta) 
real  hL-min  value  of  range  over  which  g(h)  >  0 
real  hR-  max  value  of  range  over  which  g(h)>  0 
logical  reampflag 

Computes  the  fluxes  across  each  edge  by  integrating  the  ice 
area  and  tracers  over  each  flux  triangle.  Input  variables  have 
the  same  meanings  as  in  the  main  subroutine.  Repeated  use 
of  certain  sums  makes  the  calculation  more  efficient. 

I/O:  None 
Calls:  None 

Called  by:  transport_remap 

I/O  Parameters:  real  triarea 

real  xpO,  ypO 

real  xpl,  ypl 

real  xp2,  yp2 

real  xp3,  yp3 

integer  iflux,  jflux 

real  aic,  aix,  aiy 

real  aiflx 

real  trc,  trx,  try 

real  atflx 
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Subroutine 

freeboard 


from  coupler 


frzmlt_bottom_!ateral 


get_sum 


getflux 


Description 

If  there  is  enough  snow  to  lower  the  ice/snow  interface 
below  sea  level,  this  subroutine  converts  enough  snow  to  ice 
to  bring  the  interface  back  to  sea  level. 

NOTE:  Subroutine  rebin  should  be  called  after  freeboard  to 
make  sure  ice  thicknesses  are  within  category  bounds. 

I/O:  None 
Calls:  None 
Called  by:  thennoitd 
I/O  Parameters:  None 

Reads  input  data  from  coupler  to  sea  ice  model. 

I/O:  100  (write) 

Calls:  get_sum,  ice_timer_start,  ice_timer_stop, 
cpl  interface  contractrecv,  tarea,  hm,  bound,  anglet,  t2ugrid 
Called  by:  icemodel 
I/O  Parameters:  None 

Adjusts  frzmlt  to  account  for  changes  in fhnet  since 
from  coupler.  Computes  heat  flux  to  the  bottom  surface. 
Computes  the  fraction  of  ice  that  melts  laterally. 

I/O:  None 
Calls:  None 

Called  by:  thermovertical 

I/O  Parameters:  real  Tbot-  ice  bottom  surface  temp  (°C) 
real  fbot-  heat  flux  to  ice  bottom  (W/m2) 
real  rside-  fraction  of  ice  that  melts  laterally 

Computes  a  (weighted)  sum  over  the  global  grid. 

If flag  =  1  then  workl  is  weighted  by  work2  before  being 
added  to  work3. 

I/O:  None 

Calls:  iceglobalrealsum 
Called  by:  runtime  diags,  to  coupler 
I/O  Parameters:  integer  flag 
real  workl,  work2,  work3 
real  gsum 

Gets  forcing  data  and  interpolates  as  necessary. 

I/O:  None 
Calls:  sss_sst_restore 
Called  by:  icemodel 
I/O  Parameters:  None 
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Subroutine 

global  -Conservation 


global _gather 


global _scatter 


ice  beast  char 


ice  beast  iscalar 


ice_bcast_logical 


Description 

Checks  whether  values  of  conserved  quantities  have 
changed.  An  error  probably  means  that  ghost  cells  are 
treated  incorrectly. 

I/O:  stdout 

Calls:  abort  ice 

Called  by:  transportremap 

I/O  Parameters:  real  asuminit-  initial  global  ice  area 
real  asumfinal-  final  global  ice  area 
real  atsuminit-  initial  global  ice  area*tracer 
real  atsum_final-  final  global  ice  area*  tracer 

Gathers  a  distributed  array  and  strips  off  ghost  cells  to  create 
a  local  array  with  global  dimensions. 

I/O:  None 

Calls:  icemodelsize,  iceconstants 

Called  by:  icecdf,  ice_write,  init_diags,  init_mass_diags, 

runtime  diags,  tlatlon 

I/O  Parameters:  real  work 

Scatters  a  global  array  and  adds  ghost  cells  to  create  a 
distributed  array. 

I/O:  None 

Calls:  ice  model  size,  ice  constants 

Called  by:  columngrid,  init  grid,  ice  read,  rectgrid,  tlatlon 

I/O  Parameters:  real  work 

Broadcasts  a  scalar  character  value  to  all  processors. 

I/O:  None 

Calls:  None 

Called  by:  inputdata 

I/O  Parameters:  character  charval 

Broadcasts  an  integer  scalar  character  value  to  all  processors. 
I/O:  None 
Calls:  None 

Called  by:  init_cpl,  init  hist,  restartfile 
I/O  Parameters:  integer  ival 

Broadcasts  a  scalar  logical  value  to  all  processors. 

I/O:  None 
Calls:  None 

Called  by:  ice  read,  init  hist,  input  data 
I/O  Parameters:  logical  logval 
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Subroutine 

ice  beast  rscalar 


ice_coupling_setup 


ice_global_real_minmax 


ice _g!obaI real  maxval 


ice_global_realjninval 


Description 

Broadcasts  a  real  scalar  character  value  to  all  processors. 
I/O:  None 
Calls:  None 

Called  by:  init  cpl,  input  data,  restartfile 
I/O  Parameters:  real  val 

This  routine  sets  the  model  MPI  communicators  and  task 
IDs  from  CCSM  share  code. 

I/O:  stdout 

Calls:  cpl_interface_init,  shr_sys_flush 
Called  by:  setup_mpi 

I/O  Parameters:  character  in  model  name-  input 
integer  modelcomm  -  communicator  for  model  (output) 

Detennines  and  writes  both  minimum  and  maximum  over 
global  grid  and  then  prints. 

I/O:  stdout 

Calls:  ice  fileunits,  iceglobalrealminval, 

iceglobalrealmaxval 

Called  by:  None 

I/O  Parameters:  integer  nc 

real  work(nc) 

character  string 

Computes  the  maximum  over  the  global  grid. 

I/O:  None 
Calls:  None 

Called  by:  iceglobalrealminmax,  ice_timer_print, 

runtimediags 

I/O  Parameters:  integer  nc 

real  work(nc) 

Computes  the  minimum  over  the  global  grid. 

I/O:  None 
Calls:  None 

Called  by:  ice  global  real  minmax,  ice_timer_print 
I/O  Parameters:  integer  nc 
real  work(nc) 
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Subroutine 

ice_global_reaI_sum 


ice_open 


ice  read 


Description 

Sums  a  given  array  over  the  global  grid. 

I/O:  None 
Calls:  mpiallreduce 

Called  by:  conserved_sums,  get_sum,  runtime_diags 
I/O  Parameters:  integer  nc 
real  work(nc) 

Opens  an  unformatted  file  for  reading.  The  indication  for 
whether  the  file  is  sequential  or  direct  access  is  found  in 
nbits. 

I/O:  open  nu  (unformatted  direct  access  open) 

Calls:  None 

Called  by:  dumpfile,  pipsgrid,  popgrid,  restartfile, 

read_clim_data,  read_data,  sss_clim,  sst_ic 

I/O  Parameters:  integer  nu-  unit  number 

integer  nbits-  no.  of  bits/variable  (0  for  sequential  access) 

character  filename 

Reads  an  unfonnatted  file.  Work  is  a  real  array,  atype 
indicates  the  fonnat  of  the  data. 

I/O:  stdout,  read  nu  (unfonnatted  read) 

Calls:  None 

Called  by:  pipsgrid,  popgrid,  readclimdata,  read  data, 

restartfile,  sss_clim,  sst_ic 

I/O  Parameters:  integer  nu-  unit  number 

integer  nrec-  record  number  (0  for  sequential  access) 

real  work-  output  array  (real,  8-byte) 

character  atype-  fonnat  for  input  anay  (real/int,  4-byte/8- 

byte) 

logical  scatter-  if  true,  scatter  the  data 
logical  diag-  if  true,  write  diagnostic  output 
logical  ignore  eof,  hit  eof 


98 


PIPS  3.0  SDD 


Subroutine 

ice_strength 


ice  timer  clear 


ice_timer_print 


ice  timer  start 


ice_timer_stop 


Description 

Computes  the  strength  of  the  ice  pack,  defined  as  the  energy 
(J/m2)  dissipated  per  unit  area  removed  from  the  ice  pack 
under  compression  and  assumed  proportional  to  the  change 
in  potential  energy  caused  by  ridging. 

References:  [35]  and  [16]. 

For  simpler  strength  parameterization,  see  [15]. 

I/O:  None 
Calls:  ridge_prep 
Called  by:  evp_prep 

I/O  Parameters:  integer  kstmgth-  =1  for  Rothrock 
fonnulation  [35],  0  for  Hibler  fonnation  [15] 

Initializes  timer  n  to  0.  If  n  =  -1,  all  timers  are  initialized. 

I/O:  None 
Calls:  None 
Called  by:  icemodel 

I/O  Parameters:  integer  n-  timer  number 

Prints  timing  results  of  timer  n.  If'/?  =  - 1  the  timing  results  of 
all  timers  are  printed. 

I/O:  write  100 

Calls:  icedomain,  icempiinternal,  ice  fileunits, 
iceglobalrealmaxval,  iceglobalrealminval, 

Called  by:  icemodel 

I/O  Parameters:  integer  n-  timer  number 

Begin  timing  with  timer  n. 

I/O:  None 
Calls:  timers 

Called  by:  boundijn,  evp,  icemodel,  fromcoupler, 
ridge  ice,  thennoitd,  thenno  vertical,  tocoupler, 
transportmpdata,  transportremap 
I/O  Parameters:  integer  n-  timer  number 

Ends  (or  pauses)  timing  with  timer  n. 

I/O:  None 
Calls:  timers 

Called  by:  bound  ijn,  evp,  icemodel,  from  coupler, 
ridge  ice,  thenno  itd,  thenno  vertical,  to  coupler, 
transport  mpdata,  transport  remap 
I/O  Parameters:  integer  n-  timer  number 


99 


PIPS  3.0  SDD 


Subroutine 

ice  write 


ice  write  hist 


icecdf 


init_calendar 

init_constants 

init_cpl 


Description 

Writes  an  unfonnatted  file.  Work  is  a  real  array,  atype 
indicates  the  fonnat  of  the  data. 

I/O:  stdout,  write  nu  (unformatted  write) 

Calls:  global  gather 
Called  by:  dumpfile 

I/O  Parameters:  integer  nu-  unit  number 

integer  nrec-  record  number  (0  for  sequential  access) 

real  work-  input  array 

character  atype-  format  for  output  array 

logical  gather-  if  true,  gather  the  data 

Writes  average  ice  quantities  or  snapshots. 

I/O:  None 

Calls:  icecdf,  ice_flux,  ice_albedo,  ice_mechred,  ice_grid, 

ice_calendar,  ice_state,  ice_dyn_evp,  ice_constants, 

principalstress 

Called  by:  icemodel 

I/O  Parameters:  None 

Writes  netCDF  history  file. 

I/O:  stdout,  write  ncfile 

Calls:  global_gather,  ice_model_size,  ice_constants, 
icempiintemal,  ice  grid,  ice  calendar 
Called  by:  ice_write_hist 

I/O  Parameters:  None 

Initializes  calendar  variables. 

Calls:  None 
Called  by:  icemodel 
I/O  Parameters:  None 

Initializes  constants  that  are  best  defined  at  run  time  (e.g.  pi). 
Calls:  None 
Called  by:  icemodel 
I/O  Parameters:  None 

Initializes  message  passing  between  ice  and  coupler. 

I/O:  stdout 
Calls:  None 
Called  by:  icemodel 

I/O  Parameters:  None 
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Subroutine 

init diagnostics 

initdiags 

init_evp 

init Jlux 


init Jlux_atm 


init Jlux_ocn 


Description 

Initializes  diagnostic  fields  written  to  history  files. 

I/O:  None 
Calls:  ice_state 

Called  by:  icemodel,  thermo  vertical 

I/O  Parameters:  None 

Finds  tasks  for  requested  points. 

I/O:  stdout,  2234  (write) 

Calls:  global  gather,  ice  grid,  icempiintemal 
Called  by:  icemodel 
I/O  Parameters:  None 

Initializes  parameters  needed  for  EVP  dynamics. 

I/O:  stdout  (a8,  fl2.4) 

Calls:  icecalendar,  ice  fileunits,  iceflux 
Called  by:  icemodel 
I/O  Parameters:  None 

Initializes  all  of  the  fluxes  exchanged  with  the  flux  coupler 
and  some  data  derived  fields. 

I/O:  None 

Calls:  ice  constants,  ice  flux,  initfluxatm,  init  flux  ocn 
Called  by:  icemodel 

I/O  Parameters:  None 

Initializes  all  fluxes  sent  to  the  coupler  for  use  by  the 
atmospheric  model  and  a  few  state  quantities. 

I/O:  None 

Calls:  ice_state  (only:  dice) 

Called  by:  initflux 
I/O  Parameters:  None 

Initializes  fluxes  sent  to  the  coupler  for  use  by  the  ocean 
model. 

I/O:  None 
Calls:  None 

Called  by:  init  flux,  thenno  itd 
I/O  Parameters:  None 
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Subroutine 

init_getflux 


init_grid 


init  hist 


init  itd 


init_mass_diags 

init  mechred 


Description 

Determines  the  final  year  of  the  forcing  cycle  based  on 
namelist  input.  It  initializes  the  forcing  data  filenames  and 
initializes  surface  temperature  and  salinity  from  data. 

I/O:  stdout 

Calls:  ice_history  (only:  restart),  ncar_files,  sss_clim,  sst_ic 
Called  by:  icemodel 
I/O  Parameters:  None 

Horizontal  grid  initialization  routine. 

HT{N,E}  =  cell  widths  on  {N,E}  sides  of  T  cell; 

U{LAT,LONG}  =  true  {latitude,  longitude}  of  U  points; 

D{X,Y}{T,U}  =  {x,y}  spacing  centered  at  {T,U}  points. 

I/O:  None 

Calls:  None 

Called  by:  icemodel 

I/O  Parameters:  None 

Initializes  history  files. 

I/O:  stdout,  open  nunml,  read  nunml 
Calls:  abort_ice,  ice_bcast_iscalar,  ice_bcast_logical, 
ice_constants,  ice_calendar,  ice_flux  (only:  mlt_onset, 
frz_onset),.  If  coupled,  shr_sys_mod  (only:  shr_sys  Jlush) 
Called  by:  icemodel 
I/O  Parameters:  None 

Initializes  area  fraction  and  thickness  boundaries  for  the  ITD 
model. 

I/O:  stdout 
Calls:  abort  ice 
Called  by:  icemodel 
I/O  Parameters:  None 

Computes  the  global  combined  ice  and  snow  mass  sum. 

I/O:  None 

Calls:  get_sum,  global_gather,  ice_mpi_internal,  ice_state 
Called  by:  icemodel 
I/O  Parameters:  None 

Initializes  some  variables  written  to  the  history  file. 

I/O:  None 
Calls:  None 
Called  by:  icemodel 
I/O  Parameters:  None 


102 


PIPS  3.0  SDD 


Subroutine  Description 

init  remap  Initializes  grid  quantities  used  by  remapping. 

I/O:  stdout 
Calls:  abortice,  bound 
Called  by:  icemodel 
I/O  Parameters:  None 

init_state  Initializes  state  for  the  ITD  model. 

I/O:  None 

Calls:  aggregate,  bound,  bound_aggregate,  ice_constants, 

ice_flux,  ice_grid,  ice_model_size,  ice_state, 

icethermvertical 

Called  by:  icemodel 

I/O  Parameters:  None 

init_thermo_vars  For  current  thickness  category,  this  routine  initializes  the 

thermodynamic  variables  that  are  aggregated  and  sent  to  the 
coupler,  along  with  other  fluxes  passed  among  subroutines. 
I/O:  None 
Calls:  None 

Called  by:  thermovertical 

2 

I/O  Parameters:  real  strxn-  air/ice  zonal  stress  (N/m“) 

real  stryn-  air/ice  meridional  stress  (N/m  ) 

real  Trefn-  air  temp  reference  level  (K) 

real  Qrefn-  air  speed  reference  level  (kg/kg) 

real  fsensn-  surface  downward  sensible  heat  (W/m") 

real  flatn-  surface  downward  latent  heat  (W/m") 

real  fswabsn-  SW  absorbed  by  ice  (W/nr) 

real  flwoutn-  upward  LW  at  surface  (W/m2) 

real  evapn-  flux  of  vapor,  atmosphere  to  ice  (kg/nr/s) 

real  freshn-  flux  of  water,  ice  to  ocean  (kg/m2/s) 

real  fsaltn-  flux  of  salt,  ice  to  ocean  (kg/m  /s) 

real  fhnetn-  fbot,  corrected  for  surplus  energy  (W/m2) 

real  fswthrun-  SW  through  ice  to  ocean  (W/m") 

real  fsurf-  net  flux  to  top  surface,  not  inc .fcondtop 

real  fcondtop-  downward  cond  flux  at  top  surface  (W/m  ) 

real  fcondbot-  downward  cond  flux  at  bottom  surf  (W/m2) 

real  fswint-  SW  absorbed  in  ice  interior,  below  surf  (W/m-) 

2 

real  einit-  initial  energy  of  melting  (J/m  ) 
real  efinal-  final  energy  of  melting  (J/m  ) 
real  mvap-  ice/  snow  mass  sublimated/condensed  (kg/m2) 
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Subroutine 

init  thermo  vertical 


init_vertical _profile 


input_data 


Description 

Initializes  the  vertical  profile  of  ice  salinity  and  melting 
temperature. 

I/O:  None 
Calls:  iceitd 
Called  by:  icemodel 
I/O  Parameters:  None 

Given  the  state  variables  ( vicen ,  vsnon,  eicen,  esnon,  Tsfcn ), 
this  subroutine  computes  variables  needed  for  the  vertical 
thermodynamics  ( hin ,  hsn,  qin,  qsn,  Tin,  Tsn,  Tsf). 

I/O:  stdout 

Calls:  abortice,  print  state 
Called  by:  thermovertical 

I/O  Parameters:  integer  n-  thickness  category  index 
integer  icells-  number  of  cells  with  aicen  >  puny 
integer  indexi,  indxj-  compressed  indices  for  cells  with  aicen 
>  puny 

real  hin-  ice  thickness  (m) 

real  hsn-  snow  thickness  (m) 

real  hlyr-  ice  layer  thickness 

real  hin  init-  initial  value  of  hin 

real  hsn  init-  initial  value  of  hsn 

real  qsn-  snow  enthalpy 

real  Tsn-  snow  temperature 

real  Tsf-  ice/snow  surface  temperature,  Tsfcn 

real  qin-  ice  layer  enthalpy  (J/m3) 

real  Tin-  internal  ice  layer  temperatures 

real  einit-  initial  energy  of  melting  (J/m  ) 

Namelist  variables  that  are  set  to  default  values.  They  may 
be  altered  at  run  time. 

I/O:  stdout,  read  nu  ntnl,  write  1000,  1010,  1020,  1030, 
1050,  1060,  1070,  1080,  1090,  1095 
Calls:  abort_ice,  ice  albedo,  ice  bcast  char, 
ice_bcast_iscalar,  ice_bcast_logical,  ice_bcast_rscalar, 
ice_diagnostics,  ice_mechred  (only:  kstrength,  krdgjpartic, 
krdg_redist),  ice_history,  ice_calendar,  ice_coupling  (only: 
l_coupled,  runtype ),  ice_dyn_evp,  ice_itd  (only:  kitd, 
kcatbound),  ice_ocean  (only:  ocean  mixed jce ) ,  ice_flux_in 
(only:  ycycle,  fyearjnit,  atm_data_dir,  ocn_data_dir) 
Called  by:  icemodel 
I/O  Parameters:  None 
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Subroutine 

interpcoeff 


interp_coeffjnonthly 


interpolate_data 


lateral  melt 


function  lenstr  (label) 


Description 

Computes  coefficients  for  interpolating  data  to  the  current 
time  step.  It  works  for  any  data  interval  that  divides  evenly 
into  a  365-day  year  (daily,  6-hourly,  etc.).  Uses 
intern  coef  monthly  for  monthly  data. 

I/O:  None 
Calls:  None 

Called  by:  ncarbulkdat 

I/O  Parameters:  integer  recnum-  record  number  for  current 
data  value 

integer  recslot-  spline  slot  for  current  record 
real  secant-  seconds  in  data  interval 

Computes  coefficients  for  interpolating  monthly  data  to  the 
current  time  step. 

I/O:  None 
Calls:  None 

Called  by:  ncar_bulk_dat,  sss_sst_restore 

I/O  Parameters:  integer  recslot-  slot  (1  or  2)  for  current 

record 

Linear  interpolation  subroutine. 

I/O:  None 
Calls:  None 

Called  by:  ncar_bulk_dat,  sss_sst_restore 

I/O  Parameters:  real  fielddata-  two  values  used  for 

interpolation 

real  field  -  interpolated  field 

Given  the  fraction  of  ice  melting  laterally  in  each  grid  cell 
(computed  in  subroutine  frzmlt_bottom_lateral.f),  melt  the 
ice. 

I/O:  None 

Calls:  None 

Called  by:  thennoitd 

I/O  Parameters:  real  rside-  fraction  of  ice  that  melts 
laterally 

Computes  the  length  string  by  finding  the  first  non-blank 
character  from  the  right. 

I/O:  None 
Calls:  None 

Called  by:  dumpfile,  icecdf 
I/O  Parameters:  character  label 
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Subroutine 

limited _gradient 


linear  itd 


Description 

Computes  a  limited  gradient  of  the  scalar  field  phi. 

"Limited"  means  that  we  do  not  create  new  extrema  in  phi. 
For  instance,  field  values  at  the  cell  comers  can  neither 
exceed  the  maximum  of phi(i,j)  in  the  cell  and  its  eight 
neighbors,  nor  fall  below  the  minimum. 

I/O:  None 
Calls:  None 

Called  by:  construct_fields 

I/O  Parameters:  real  phi-  input  tracer  field  (mean  values  in 
each  grid  cell) 

real  cnx-  x-coordinate  of  phi  relative  to  geometric  center  of 
cell 

real  cny-  y-coordinate  of  phi  relative  to  geometric  center  of 
cell 

real  phimask-  phimask(ij)  =  1  if phi(i,j)  has  physical 
meaning,  =  0  otherwise.  For  instance,  aice  has  no  physical 
meaning  on  land  points,  and  hice  no  physical  meaning  where 
aice  =  0 

real  gx-  limited  x-direction  gradient 
real  gy-  limited  y-direction  gradient 

Ice  thickness  distribution  scheme  that  shifts  ice  among 
categories.  The  default  scheme  is  linear  remapping,  which 
works  as  follows: 

Using  the  thennodynamic  "velocities",  it  interpolates  to  find 
the  velocities  in  thickness  space  at  the  category  boundaries 
and  computes  the  new  locations  of  the  boundaries.  Then  for 
each  category,  the  scheme  computes  the  thickness 
distribution  function,  g(h),  between  hL  and  hR,  the  left  and 
right  boundaries  of  the  category. 

I/O:  stdout 

Calls:  aggregate_area,  column_sum, 
columnconservationcheck,  fit  line,  shift  ice 
Called  by:  thennoitd 

I/O  Parameters:  real  hicen  old-  starting  value  of  hicen 
real  hicen-  ice  thickness  for  each  category  (m) 
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Subroutine  Description 

load  tracers  Loads  tracer  array  and  computes  tracer  dependency  vectors. 

The  default  version  assumes  that  the  advected  tracers  are 
hice,  hsno,  Tsfc,  qice(l:nilyr),  and  qsno.  This  subroutine 
must  be  modified  if  a  different  set  of  tracers  is  to  be 
transported.  The  rule  for  ordering  tracers  is  that  a  dependent 
tracer  (such  as  qice )  must  have  a  larger  tracer  index  than  the 
tracer  it  depends  on  (i.e.,  hice). 

I/O:  None 
Calls:  None 

Called  by:  transportremap 

I/O  Parameters:  integer  n-  ice  category  index 

real  tnn-  mean  tracer  values  in  each  grid  cell 

local  maxjnin  At  each  grid  point,  this  subroutine  computes  the  local  max 

and  min  of  a  scalar  field  phi :  i.e.,  the  max  and  min  values  in 
the  nine-cell  region  consisting  of  the  home  cell  and  its  eight 
neighbors,  plus  the  neighbors  of  the  neighbors  (25  cells  in 
all). 

I/O:  None 
Calls:  bound 

Called  by:  transport_remap 
I/O  Parameters:  real  aimask 
real  tnn 
real  trmask 

real  tmin-  local  min  tracer 
real  tmax-  local  max  tracer 
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Subroutine 

locate  triangles 


make  masks 


makemask 


Description 

Computes  areas  and  vertices  of  flux  triangles  for  east  and 
north  cell  edges. 

I/O:  None 
Calls:  None 

Called  by:  transportremap 

I/O  Parameters:  real  dpx,dpy-  x,  y  coordinates  of  departure 
points  at  cell  comers 

real  triarea_e-  area  of  east-edge  flux  triangle 

triarean-  area  of  north-edge  flux  triangle 

real  xpO  e,  ypO  e-  coordinates  of  special  triangle  points 

real  xp0_e,  yp0_e-  e  for  east  edges,  n  for  north  edges 

real  xpl_e,  ypl_e 

real  xp2_e,  yp2_e 

real  xp3_e,  yp3_e 

real  xpOn,  ypOn 

real  xpln,  ypl  n 

real  xp2_n,  yp2_n 

real  xp3_n,  yp3_n 

Creates  area  and  tracer  masks.  If  an  area  is  masked  out  ( aim 
<  puny),  then  the  values  of  tracers  in  that  grid  cell  are 
assumed  to  have  no  physical  meaning.  Similarly,  if  a  tracer 
with  dependents  is  masked  out  ( abs(trm )  <  puny),  then  the 
values  of  its  dependent  tracers  in  that  grid  cell  are  assumed 
to  have  no  physical  meaning.  For  example,  the  enthalpy 
value  has  no  meaning  if  the  thickness  is  zero. 

I/O:  None 
Calls:  None 

Called  by:  transport  remap 

I/O  Parameters:  real  aim-  mean  ice  area  in  each  grid  cell 
real  aimask-  1 .  if  ice  is  present,  otherwise  =  0 
real  tnn-  mean  tracer  values  in  each  grid  cell 
real  trmask-  1 .  if  tracer  is  present,  else  =  0 

Sets  the  boundary  values  for  the  T  cell  land  mask  {km)  and 
makes  the  logical  land  masks  for  T  and  U  cells  ( tmask , 
umask).  This  routine  also  creates  hemisphere  masks  ( mask-n 
Northern,  mask-s  Southern). 

I/O:  None 
Calls:  bound 
Called  by:  init_grid 
I/O  Parameters:  None 
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Subroutine 

merge Jluxes 


mixed Jayer 


mpdata  (narrays.phi) 


Description 

Aggregates  flux  information  from  all  ice  thickness 
categories. 

I/O:  None 

Calls:  ice_state 

Called  by:  thermovertical 

I/O  Parameters:  integer  n-  thickness  category  index 

real  strxn-  air/ice  zonal  stress  (N/nr) 

real  stryn-  air/ice  meridional  stress  (N/nr) 

real  fsensn-  sensible  heat  flux  (W/nT) 

real  flatn-  latent  heat  flux  (W/m“) 

2 

real  fswabsn-  shortwave  absorbed  heat  flux  (W/m“) 

real  flwoutn-  upward  low  emitted  heat  flux  (W/m“) 

real  evapn-  evaporation  (kg/m  /s) 

real  Trefn-  air  temp  reference  level  (K) 

real  Qrefn-  air  speed  humidity  reference  level  (kg/kg) 

real  freshn-  fresh  water  flux  to  ocean  (kg/m2/s) 

real  fsaltn-  salt  flux  to  ocean  (kg/nr/s) 

real  fhnetn-  actual  ocean/ice  heat  flux  (W/m2) 

real  fswthrun-  SW  radiation  through  ice  bottom  (W/m  ) 

Computes  the  mixed  layer  heat  balance  and  updates  the  SST. 
This  routine  also  computes  the  energy  available  to  freeze  or 
melt  ice. 

NOTE:  SST  changes  due  to  fluxes  through  the  ice  are 
computed  in  ice  therm  vertical. 

I/O:  None 

Calls:  atmoboundarylayer,  ice  flux,  icecalendar  (only: 
dt),  ice_grid  (only:  Mask),  ice_state,  ice_albedo 
Called  by:  icemodel 
I/O  Parameters:  None 

Advection  according  to  mpdata. 

Reference:  [39]. 

I/O:  stdout 

Calls:  bound,  boundnarr,  ice  calendar,  ice  state  (only: 

uvel,  wel),  abort_ice 

Called  by:  transportmpdata 

I/O  Parameters:  integer  narrays 

real  phi 
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Subroutine 

NCAR  bulk  dat 


NCAR  Jiles 


pipsgrid 


popgrid 


Description 

Reads  NCAR  bulk  atmospheric  data. 

I/O:  readm,  read6,  read 
Calls:  interpcoeffmonthly,  readclimdata 
Called  by:  getflux 
I/O  Parameters:  None 

This  subroutine  is  based  on  the  LANL  file  naming 
conventions.  The  user  must  edit  it  for  other  directory 
structures  or  filenames. 

NOTE:  The  year  number  in  these  filenames  does  not  matter, 
because  subroutine  fdejyear  will  insert  the  correct  year. 

I/O:  stdout 
Calls:  file_year 
Called  by:  init_getflux 

I/O  Parameters:  integer  yr-  current  forcing  year 
PIPS  3.0  rotated  spherical  grid  and  land  mask. 


Rec  No. 

Field 

Units 

Land  Mask  1 

KMT 

1  grid  file  with 
real  KMT's  grid 

2 

ULAT 

radians 

3 

ULON 

radians 

4 

HTN 

cm 

5 

HTE 

cm 

6 

HUS 

cm 

7 

HUW 

cm 

8 

ANGLE 

radians 

I/O:  stdout 
Calls:  ice_read_write 
Called  by:  init_grid 
I/O  Parameters:  None 

Reads  and  sets  POP  displaced  pole  grid  and  land  mask. 
See  Subroutine  pipsgrid  description  above  for  record 
number,  field  and  units,  as  they  are  identical. 

I/O:  None 

Calls:  ice_read_write,  ice_open,  ice_read 
Called  by:  init_grid 
I/O  Parameters:  None 
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Subroutine 

prepare Jorcing 

principal _stress 


print_state 


read  clim  data 


Description 

Finishes  the  task  of  manipulating  forcing. 

I/O:  None 

Calls:  ice  state  (only:  aice )-  lipscombtune 
Called  by:  getflux 
I/O  Parameters:  None 

Computes  principal  stresses  for  comparison  with  the 
theoretical  yield  curve;  northeast  values. 

I/O:  None 
Calls:  None 

Called  by:  ice_write_hist 

I/O  Parameters:  None 

Prints  ice  state  for  a  specified  grid  point.  This  routine  is 
useful  for  debugging. 

I/O:  stdout 

Calls:  icemodelsize,  icekindsmod,  ice  state,  iceitd, 
iceflux 

Called  by:  conservationcheckvthermo,  debug  ice, 
init_vertical_profile,  temperaturechanges 
I/O  Parameters:  character  plabel-  input 
character  i,j-  input 

Reads  annual  climatological  data  needed  for  interpolation,  as 
in  read_data.  It  assumes  a  one-year  cycle  of  climatological 
data,  so  that  there  is  no  need  to  get  data  from  other  years  or 
to  extrapolate  data  beyond  the  forcing  time  period. 

I/O:  stdout 

Calls:  ice_diagnostics  (for  debugging),  ice_open,  ice_read 

Called  by:  ncar_bulk_dat,  sss_sst_restore 

I/O  Parameters:  logical  readflag 

integer  reed-  baseline  record  number 

integer  imx,  ixx,  ipx-  record  numbers  of  three  data  values 

relative  to  reed 

character  data  file 

real  field  data-  two  values  needed  for  interpolation 
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Subroutine 

read  data 


rebin 


rectgrid 


Description 

Reads  data  for  interpolation.  If  data  is  at  the  beginning  of  a 
one-year  record,  this  subroutine  gets  data  from  the  previous 
year.  If  data  is  at  the  end  of  a  one-year  record,  it  gets  data 
from  the  following  year.  If  no  earlier  data  exists  (beginning 
of  Jyearinit),  then: 

(1)  For  monthly  data,  get  data  from  the  end  of  yearjinal. 

(2)  For  more  frequent  data,  let  the  imx  value  equal  the 
first  value  of  the  year. 

If  no  later  data  exists  (end  of JyearJinal),  then: 

(1)  For  monthly  data,  get  data  from  the  beginning  of 
Jyearinit. 

(2)  For  more  frequent  data,  let  the  ipx  value  equal  the  last 
value  of  the  year. 

In  other  words,  assume  persistence  when  daily  or  6-hourly 
data  is  missing,  and  assume  periodicity  when  monthly  data  is 
missing. 

I/O:  stdout 

Calls:  file_year,  ice_diagnostics,  ice_open,  ice_read 

Called  by:  ncarbulkdat 

I/O  Parameters:  logical  flag 

integer  reed-  baseline  record  number 

integer  yr-  year  of  forcing  data 

integer  imx,  ixx,  ipx-  record  numbers  of  three  data  values 
relative  to  reed 

integer  maxrec-  maximum  record  value 

real  field  data-  two  values  needed  for  interpolation 

Rebins  thicknesses  into  defined  categories. 

I/O:  None 

Calls:  ice_grid,  shift_ice 
Called  by:  thennoitd 
I/O  Parameters:  None 

Regular  rectangular  grid  and  mask. 

I/O:  None 

Calls:  global_scatter,  ice_model_size 
Called  by:  init_grid 
I/O  Parameters:  None 
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Subroutine 

reduce  area 


restartfile 


ridge Jce 


Description 

Reduces  area  when  ice  melts  for  special  case  ncat=l.  Use 
a  CSM  1 .0-like  method  of  reducing  ice  area  when  melting 
occurs.  Assumes  only  half  the  ice  volume  change  goes  to 
thickness  decrease  and  the  other  half  to  reduction  in  ice 
fraction. 

I/O:  None 
Calls:  icegrid 
Called  by:  thennoitd 

I/O  Parameters:  real  hicelold-  old  ice  thickness  for 
category  1  (m) 

real  hicel-  new  ice  thickness  for  category  1  (in) 

Restarts  from  a  dumpfile. 

I/O:  read,  write  nu_rst_pointer,  write  filename,  read 
nu  restart,  stdout 

Calls:  aggregate,  bound,  bound_aggregate,  bound_state, 
ice_bcast_iscalar,  ice_bcast_rscalar,  ice_model_size, 
ice_flux,  ice_mpi_intemal,  ice_grid,  ice_calendar,  ice_state, 
ice_dyn_evp,  ice_itd,  ice_ocean,  ice_open,  ice_read,  (only: 
ocean  mixed  jce ) ,  ice_coupling  (only:  l_coupled),  lenstr 
Called  by:  icemodel 
I/O  Parameters:  None 

Computes  changes  in  the  ice  thickness  distribution  due  to 
divergence  and  shear. 

References:  [13],  [16],  [35],  [43]. 

I/O:  stdout 

Calls:  abort_ice,  asum_ridging,  ice_timers,  ice_timer_start, 
ice_timer_stop,  ice_flux,  ridge_prep,  ridge_shift 
Called  by:  icemodel 

I/O  Parameters:  real  Delta-  in  the  denominator  of  zeta,  eta 

(1/s) 

real  divu-  strain  rate  I  component,  velocity  divergence  (1/s) 
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Subroutine 

ridge _prep 


ridge_shift 


runtime  diags 


scale Jluxes 


Description 

Preparation  for  ridging  and  strength  calculations.  Computes 
the  thickness  distribution  of  the  ice  and  open  water 
participating  in  ridging  and  of  the  resulting  ridges. 

This  version  includes  new  options  for  ridging  participation 

and  redistribution.  The  new  participation  scheme  improves 

model  stability  by  increasing  the  time  scale  for  large  changes 

in  ice  strength.  The  new  redistribution  scheme  improves  the 

agreement  between  ITDs  of  modeled  and  observed  ridges. 

I/O:  None 

Calls:  asumridging 

Called  by:  ice_strength,  ridge_ice 

I/O  Parameters:  None 

Shifts  ridging  ice  among  thickness  categories.  It  removes 
area,  volume,  and  energy  from  each  ridging  category  and 
adds  them  to  thicker  ice  categories. 

I/O:  stdout 

Calls:  abort  ice,  columnconservationcheck,  columnsum 
Called  by:  ridge  ice 

I/O  Parameters:  real  opening-  rate  of  opening  due  to 
divergence/sheer 

real  closing  gross-  rate  at  which  area  is  removed,  not 
counting  area  of  new  ridges 

real  msnow_mlt-  mass  of  snow  added  to  ocean  (kg/m“) 
real  esnow_mlt-  energy  needed  to  melt  snow  in  ocean  (J/m2) 

Writes  diagnostic  infonnation,  such  as  max,  min,  global 
sums,  etc.,  to  standard  out. 

I/O:  stdout,  799,  800,  801,  899,  900,  901,  902,  903  (write) 
Calls:  global_gather,  ice_global_real_maxval, 
ice  model  size,  ice  flux,  icealbedo,  iceglobalrealsum, 
icempiintemal,  ice  state,  ice  itd,  get  sum.  If  coupled, 
shr_sys_mod  (only:  shr_sys  Jlush ) 

Called  by:  icemodel 
I/O  Parameters:  None 

Divides  ice  fluxes  by  ice  area  before  sending  them  to  the 
coupler,  since  the  coupler  multiplies  by  ice  area.  This  is  the 
ice  area  at  the  beginning  of  the  timestep,  i.e.  the  value  sent  to 
the  coupler. 

I/O:  None 
Calls:  icealbedo 
Called  by:  icemodel 
I/O  Parameters:  None 
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Subroutine 

scale  hist Jluxes 


setup jnpi 


shift Jce 


sss  dim 


sss  sst  restore 


Description 

Divides  ice  fluxes  by  ice  area  used  by  the  coupler  before 
writing  out  diagnostics.  Aice  init  is  the  ice  area  saved  from 
coupling.  This  makes  the  fluxes  written  to  the  history  file 
consistent  with  those  sent  to  the  coupler. 

I/O:  None 
Calls:  None 
Called  by:  icemodel 
I/O  Parameters:  None 

This  routine  initializes  MPI  for  either  internal  parallel 
processing  or  for  message  passing  with  the  coupler. 

I/O:  stdout 

Calls:  abortice,  icempiintemal,  icecoupling, 

icecouplingsetup 

Called  by:  icemodel 

I/O  Parameters:  None 

Shifts  ice  across  category  boundaries,  conserving  area, 
volume,  and  energy. 

I/O:  stdout 

Calls:  abort_ice,  ice_flux,  ice_work  (only:  worka ) 

Called  by:  linear  itd,  rebin 

I/O  Parameters:  integer  donor-  donor  category  index 
real  daice-  ice  area  transferred  across  boundary 
real  dvice-  ice  volume  transferred  across  boundary 
real  hicen-  ice  thickness  for  each  category  (m) 

NOTE:  Third  index  of  donor ;  daice,  dvice  should  be  ncat-1, 
except  that  compilers  would  have  trouble  when  neat  =  1. 

Creates  an  annual  mean  climatology  for  Levitus  SSS  from  a 
12-month  climatology. 

I/O:  stdout 

Calls:  ice_open,  ice_read,  ice_work  (only:  worka) 

Called  by:  init_getflux 
I/O  Parameters:  None 

Interpolates  monthly  SSS  and  SST  data  to  timestep.  This 
subroutine  restores  SST  computed  by  the  ice  model  to  data. 
I/O:  stdout,  readm 

Calls:  completegetfluxocn,  interpcoeffmonthly, 
interpolate  data,  readclimdata 
Called  by:  getflux 
I/O  Parameters:  None 
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Subroutine 

sst  ic 


stepu 


stress 


surface Jluxes 


Description 

Reads  SST  data  for  current  month,  and  adjusts  SST  based  on 
freezing  temperature.  This  routine  does  not  interpolate. 

I/O:  stdout 

Calls:  iceopen,  ice  read 
Called  by:  init_getflux 
I/O  Parameters:  None 

Calculation  of  the  surface  stresses  and  integration  of  the 
momentum  equation  to  find  velocity  (u,v). 

I/O:  None 
Calls:  iceflux 
Called  by:  evp 
I/O  Parameters:  None 

Computes  the  rates  of  strain  and  internal  stress  components 
for  each  of  the  four  comers  on  each  T-grid  cell. 

I/O:  None 
Calls:  None 
Called  by:  evp 

I/O  Parameters:  integer  ksub-  subcycling  step  input) 

Computes  radiative  and  turbulent  fluxes  and  their  derivatives 
with  respect  to  Tsf 
I/O:  None 
Calls:  None 

Called  by:  temperature_changes 

I/O  Parameters:  integer  isolve-  no.  of  cells  with  temps  not 
converged 

integer  indxii,  indxjj-  compressed  indices  for  cells  not 
converged 

real  Tsf-  ice/snow  surface  temperature,  Tsfcn 

real  fswsfc-  SW  absorbed  at  ice/  snow  surface  (W/m2) 

real  fsensn-  surface  downward  sensible  heat  (W/m  ) 

real  flatn-  surface  downward  latent  heat  (W/m") 

real  flwoutn-  upward  LW  at  surface  (W/m2) 

real  fsurf-  net  flux  to  top  surface,  not  incl.  fcondtop 

real  dfsens_dT-  derivative  of fsens  with  respect  to  7)/  (W/m 2 

/deg) 

real  dflat_dT-  deriv  of flat  with  respect  to  ?V(W/m2  /deg) 
real  dflwout_dT-  deriv  of flwout  with  respect  to  Ts/( W/m2 
/deg) 

real  dfsurf_dT-  derivative  of fsurf  with  respect  to  Tsf 
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Subroutine 

t2ugrid 


temperature  changes 


Description 

Transfers  from  T-cell  centers  to  U-cell  centers.  Writes  work 
into  another  array  that  has  ghost  cells. 

I/O:  None 

Calls:  bound,  to  ugrid 

Called  by:  evp_prep,  fromcoupler 

I/O  Parameters:  real  work 

Computes  new  surface  temperature  and  internal  ice  and 
snow  temperatures.  It  includes  effects  of  salinity  on  sea  ice 
heat  capacity  in  a  way  that  conserves  energy  [7].  New 
temperatures  are  computed  iteratively  by  solving  a 
tridiagonal  system  of  equations.  Heat  capacity  is  updated 
with  each  iteration.  Finite  differencing  is  backward  implicit. 
I/O:  stdout 

Calls:  abortice,  absorbedsolar,  conductivity,  print  state, 
surface_fluxes,  tridiag_solver 
Called  by:  thennovertical 

I/O  Parameters:  integer  n-  thickness  category  index 
integer  icells-  number  of  cells  with  aicen  >  puny 
integer  indxi,  indxj-  compressed  indices  for  cells  with  aicen 
>  puny 

real  hlyr-  ice  layer  thickness 

real  hsn-  snow  thickness  (in) 

real  Tbot-  ice  bottom  surface  temp  (°  C) 

real  fbot-  ice-ocean  heat  flux  at  bottom  surface  (W/m2) 

real  qsn-  snow  enthalpy 

real  Tsn-  internal  snow  temperature 

real  Tsf-  ice/snow  surface  temperature,  Tsfcn 

real  fsensn-  surface  downward  sensible  heat  (W/m2) 

real  flatn-  surface  downward  latent  heat  (W/m“) 

real  fswabsn-  shortwave  absorbed  by  ice  (W/nr) 

real  flwoutn-  upward  LW  at  surface  (W/m') 

real  fswthrun-  SW  through  ice  to  ocean  (W/m2) 

real  fhnetn-  fbot,  corrected  for  any  surplus  energy 

real  fsurf-  net  flux  to  top  surface,  not  including  fcondtop 

real  fcondtop-  downward  cond  flux  at  top  surface  (W/m  ) 

real  fcondbot-  downward  cond  flux  at  bottom  surface 

(W/m2) 

real  fswint-  SW  absorbed  in  ice  interior,  below  surface 
(W/m2) 

real  qin-  ice  layer  enthalpy  (J/m  ) 
real  Tin-  internal  ice  layer  temperatures 
real  einit-  initial  energy  of  melting  (J/m2) 


117 


PIPS  3.0  SDD 


Subroutine  Description 

thermo _itd  Driver  for  post-coupler  thermodynamic  changes  not  needed 

for  coupling:  transport  in  thickness  space,  lateral  growth  and 
melting,  and  freeboard  adjustment. 

NOTE:  Ocean  fluxes  are  initialized  here. 

I/O:  None 

Calls:  add_new_ice,  aggregate,  freeboard,  ice_timers, 
ice_itd_linear,  ice_thenn_vertical,  ice_timer_start, 
icetimerstop,  init  flux  ocn,  lateralmelt,  linear  itd, 
reduce  area,  rebin,  zapsmallareas 
Called  by:  icemodel 
I/O  Parameters:  None 

thermo  vertical  Driver  for  updating  ice  and  snow  internal  temperatures  and 

computing  thennodynamic  growth  rates  and  atmospheric 
fluxes. 

I/O:  None 

Calls:  add_new_snow,  atmo_boundary_layer, 
conservation  check  vthenno,  frzmltbottomlateral, 
ice_atmo,  ice_timers,  ice_timer_start,  ice_timer_stop, 
ice_ocean,  ice_work  (only:  worka,  workb ),  init_diagnostics, 
initfluxatm,  init  thenno  vars,  init_vertical_profile, 
merge_fluxes,  temperature_changes,  thickness_changes, 
updatestatevthenno 
Called  by:  icemodel 
I/O  Parameters:  None 
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Subroutine 

thickness  changes 


timers 


tlatlon 


tocoupler 


Description 

Computes  growth  and/or  melting  at  the  top  and  bottom 
surfaces. 

I/O:  None 
Calls:  None 

Called  by:  thennovertical 

I/O  Parameters:  integer  n-  thickness  category  index 
integer  icells-  number  of  cells  with  aicen  >  puny 
integer  indxi,  indxj-  compressed  indices  for  cells  with  aicen 
>  puny 

real  fbot  -  ice-ocean  heat  flux  at  bottom  surface  (W/m2) 

real  Tbot  -  ice  bottom  surface  temperature  (°  C) 

real  flatn  -  surface  downward  latent  heat  (W/m2) 

real  fsurf  -  net  flux  to  top  surface,  not  including  fcondtop 

real  fcondtop-  downward  cond  flux  at  top  surface  (W/m2) 

real  fcondbot-  downward  cond  flux  at  bottom  surface 

(W/m2) 

real  qin-  ice  layer  enthalpy  (J/m  ) 

real  fhnetn-  fbot,  corrected  for  any  surplus  energy 

real  hlyr  -  ice  layer  thickness 

real  hsn  -  snow  thickness  (m) 

real  qsn  -  snow  enthalpy 

real  efinal-  final  energy  of  melting  (J/m2) 

real  mvap  -  ice/snow  mass  sublimated/condensed  (kg/nr) 

real  hin-  total  ice  thickness 

Does  the  work. 

I/O:  None 
Calls:  None 

Called  by:  ice  timer  start,  icetimerstop 
I/O  Parameters:  real  tl 

Initializes  latitude  and  longitude  on  T  grid. 

I/O:  stdout 

Calls:  bound,  global  gather,  global  scatter,  ice  model  size, 
ice_read_write  (if  reading  ULAT,  ULON  directly  from  file) 
Called  by:  init_grid 
I/O  Parameters:  None 

Sends  data  from  PIPS  3.0  model  to  coupler. 

I/O:  stdout,  100  (write) 

Calls:  get_sum,  ice_timer_start,  ice_timer_stop 
Called  by:  initcpl 
I/O  Parameters:  None 
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Subroutine 

tojgrid 


tougrid 


transport jnpdata 


transport remap 


Description 

Shifts  quantities  from  the  U-cell  midpoint  (work/)  to  the  T- 
cell  midpoint  (work!). 

I/O:  None 
Calls:  None 
Called  by:  u2tgrid 

I/O  Parameters:  real  workl,  real  work2 

Shifts  quantities  from  the  T-cell  midpoint  (workl)  to  the  U- 
cell  midpoint  (work!). 

I/O:  None 
Calls:  None 

Called  by:  evp_prep,  t2ugrid 

I/O  Parameters:  real  workl,  real  work2 

Computes  the  transport  equations  for  one  timestep  using 
mpdata.  Sets  several  fields  into  a  work  array  and  passes  it  to 
mpdata  routine. 

I/O:  stdout 

Calls:  boundnarr,  check  state,  ice  flux,  ice  timers, 
ice_state,  ice_itd  (check_state),  ice_timer_start, 
icetimerstop,  mpdata 
Called  by:  icemodel 
I/O  Parameters:  None 

This  subroutine  solves  the  transport  equations  for  one 
timestep  using  the  conservative  remapping  scheme 
developed  by  John  Dukowicz  and  John  Baumgardner  (DB) 
and  modified  for  sea  ice  by  William  Lipscomb  and  Elizabeth 
Hunke.  This  scheme  is  second-order  accurate,  except  where 
gradients  are  limited  to  preserve  monotonicity.  It  is 
compatible  for  tracers;  that  is,  it  does  not  produce  new 
extrema  in  thickness  or  enthalpy. 

References:  [9],  [26]. 

I/O:  None 

Calls:  bound,  bound  narr,  bound  state,  bound  sw, 
checkmonotonicity,  conserved  sums,  construct_fields, 
departure_points,  flux  integrals,  global  conservation, 
ice_timer_start,  ice_timer_stop,  load_tracers, 
localmaxmin,  locate  triangles,  make  masks, 
unload  tracers,  update  fields 
Called  by:  icemodel 
I/O  Parameters:  None 
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Subroutine 

triangle -Coordinates 


tridiag_solver 


Description 

For  each  triangle,  this  subroutine  finds  the  coordinates  of  the 
four  points  needed  to  compute  integrals  of  cubic 
polynomials,  using  a  fonnula  from  [35].  (Section  8.8, 
fonnula  3.1.)  Quadratic  functions  can  be  integrated  using  3- 
point  fonnulas,  but  it  is  more  efficient  to  use  a  single 
fonnula  for  both  quadratics  and  cubics. 

The  fonnula  is  as  follows: 

13  =  integral  of  f [x,y)*dA 
=  AR*  [-9/16  *  f(xO,  yO) 

+  25/48  *  (f \xl,yl)  +  f(x2,y2)  +  f(x3,y3))] 
where  13  is  the  integral  of  a  polynomial  of  3rd  order  or 
lower,  AR  is  the  area  of  the  triangle,  (x0,y0)  is  the  midpoint, 
and  the  other  three  points  are  located  2/5  of  the  way  from 
the  midpoint  to  each  of  the  three  vertices. 

I/O:  None 
Calls:  None 

Called  by:  transportremap 
I/O  Parameters:  real  triarea-  areas  of  flux  triangles 
real  xpO,  ypO-  coordinates  of  triangle  points 
xp 1 ,  yp 1 

xp2, yp2 
xp3,  yp3 

Tridiagonal  matrix  solver.  It  is  used  to  solve  the  implicit 
vertical  heat  equation  in  ice  and  snow. 

I/O:  None 
Calls:  None 

Called  by:  temperature_changes 

I/O  Parameters:  integer  isolve-  number  of  cells  with  temps 
not  converged 

integer  indxii,  indxjj-  compressed  indices  for  cells  not 
converged 

integer  nmat-  matrix  dimension 
real  diag-  diagonal  matrix  elements 
real  sbdiag-  sub-diagonal  matrix  elements 
real  spdiag-  super-diagonal  matrix  elements 
real  rhs-  rhs  of  tri-diagonal  matrix  equation 
real  xout-  solution  vector 
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Subroutine 

u2tgrid 


unload  tracers 


update Jields 


Description 

Transfers  from  U-cell  centers  to  T-cell  centers.  Writes  work 
into  another  array  that  has  ghost  cells. 

I/O:  None 

Calls:  bound,  to  tgrid 
Called  by:  evpfinish 
I/O  Parameters:  real  work 

Converts  from  tracer  array  back  to  state  variables.  Like 
loadjracers,  this  subroutine  must  be  modified  if  a  different 
set  of  tracers  is  to  be  transported. 

I/O:  None 
Calls:  None 

Called  by:  transportremap 

I/O  Parameters:  integer  n-  ice  category  index 

real  tnn-  mean  tracer  values  in  each  grid  cell 

Given  fluxes  through  cell  edges,  this  subroutine  computes 
new  area  and  tracers. 

I/O:  None 

Calls:  abort  ice 

Called  by:  transport  remap 

I/O  Parameters:  real  aiflxe,  aiflxn-  flux  of  a  ice  through  east 
and  north  cell  edges 
real  aim-  mean  ice  area 

real  atflxe,  atflxn-  flux  of  r//ce*  tracer  through  E  and  N  cell 
edges 

real  tnn-  mean  tracers 

integer  n-  ice  category  index  (for  error  diagnostics) 
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Subroutine 

update_state_vthermo 


zap  small  areas 


Description 

Given  the  vertical  thenno  state  variables  (hin,  hsn,  qin, 
qsn,  Tsf),  this  subroutine  computes  the  new  ice  state 
variables  ( vicen ,  vsnon,  eicen,  esnon,  Tsfcii). 

State  variables  are  zeroed  out  if  ice  has  melted  entirely. 

I/O:  None 
Calls:  None 

Called  by:  thermovertical 

I/O  Parameters:  integer  n-  thickness  category  index 
integer  icells-  number  of  cells  with  aicen  >  puny 
integer  indxi,  indxj-  compressed  indices  for  cells  with  aicen 
>  puny 

real  hin-  ice  thickness  (m) 

real  hsn-  snow  thickness  (m) 

real  qsn-  snow  enthalpy 

real  Tsf-  ice/snow  surface  temperature,  Ts/cn 

real  qin-  ice  layer  enthalpy  (J/m3) 

Eliminates  very  small  ice  areas. 

I/O:  stdout 

Calls:  abort_ice,  ice_flux,  ice_calendar  (only:  dt) 

Called  by:  icemodel,  thenno  itd 

I/O  Parameters:  None 
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6.0  PIPS  3.0  Primary  Variables  and  Parameters 

The  following  table  defines  many  of  the  symbols  frequently  used  in  the  PIPS  3.0  code.  Values 
appearing  in  this  list  are  either  fixed  or  recommended;  most  namelist  parameters  are  indicated  (*) 
with  their  default  values.  For  other  namelist  options,  see  Table  2.  All  quantities  in  the  code  are 
expressed  in  MKS  units  (temperatures  may  take  either  Celsius  or  Kelvin  units). 


Name 

Description 

Default  Values 

A 

advection 

type  of  advection  algorithm  used 

‘remap’ 

ahmax 

thickness  above  which  ice  albedo  is  constant 

0.5  m 

aiceO 

fractional  open  water  area 

aice(n) 

total  concentration  of  ice  in  grid  cell  (in  category  n) 

aice  init 

concentration  of  ice  at  beginning  of  dt  (for 
diagnostics) 

ain  min 

minimum  fractional  ice  area  allowed  in  each 
category 

albicei 

*near  infrared  ice  albedo  for  thicker  ice 

0.36 

albicev 

*visible  ice  albedo  for  thicker  ice 

0.78 

albsnowi 

*near  infrared,  cold  snow  albedo 

0.70 

albsnowv 

*visible,  cold  snow  albedo 

0.98 

albocn 

ocean  albedo 

0.06 

alpha 

floe  shape  constant  for  lateral  melt 

0.66 

astar 

e-folding  scale  for  participation  function 

0.05 

awtidf 

weighting  factor  for  near-ir,  diffuse  albedo 

0.16 

awtidr 

weighting  factor  for  near-ir,  direct  albedo 

0.31 

awtvdf 

weighting  factor  for  visible,  diffuse  albedo 

0.24 

awtvdr 

weighing  factor  for  visible,  direct  albedo 

0.29 

ANGLE 

for  conversions  between  the  ocean  grid  and  lat-lon 
grids 

ANGLET 

ANGLE  converted  to  T-cells 

atm  data  dir 

*directory  for  atmospheric  forcing  data 

avgsiz 

number  of  fields  that  may  be  written  to  history  file 

91 

c 
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Name 

Description 

Default  Values 

Cf 

ratio  of  ridging  work  to  PE  change  in  ridging 

17. 

char  len 

length  of  character  variable  strings 

80 

char  len  long 

length  of  longer  character  variable  strings 

128 

check_step 

time  step  for  writing  debugging  data 

cldf 

cloud  fraction 

congel 

basal  ice  growth 

m 

cosw 

cosine  of  the  turning  angle  in  water 

1. 

cpair 

specific  heat  of  air 

1005.0  J/kg/K 

cp_wv 

specific  heat  of  water  vapor 

1.81  x  103  J/kg/K 

cp_ice 

specific  heat  of  fresh  ice 

2106.  J/kg/K 

cpocn 

specific  heat  of  sea  water 

4218.  J/kg/K 

cm  to  m 

cm  to  meters  conversion 

0.01 

c  <n> 

rca](/?) 

Cs 

fraction  of  shear  energy  contributing  to  ridging 

0.5 

Cstar 

constant  in  Hibler  ice  strength  formula 

20 

D 

daidtd 

ice  area  tendency  due  to  dynamics/transport 

1/s 

daidtt 

ice  area  tendency  due  to  thermodynamics 

1/s 

dalb  mlt 

[see  icealbedo.F] 

-0.075 

dalb  mlti 

[see  ice  albedo.F] 

-0.100 

dalb  mltv 

[see  ice  albedo.F] 

-0.150 

dardgldt 

rate  of  fractional  area  loss  by  ridging  ice 

1/s 

dardg2dt 

rate  of  fractional  area  gain  by  new  ridges 

1/s 

dvirdgdt 

ice  volume  ridging  rate 

m/s 

dblkind 

definition  of  double  precision 

selected  real  kind 
(13) 

dbug 

*write  forcing  data  diagnostics 

.false. 

Delta 

function  of  strain  rates  (see  Section  5. 2.2. 4) 

depressT 

ratio  of  freezing  temps  to  salinity  of  brine 

0.054  deg/psu 

diagfile 

*diagnostic  output  file  (alternative  to  stdout) 

diag_type 

*where  diagnostic  output  is  written 

stdout 
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Name 

Description 

Default  Values 

diagfreq 

*how  often  diagnostic  output  is  written  (10  = 
once/10  dt 

24 

divu 

strain  rate  I  component,  velocity  divergence 

1/s 

divuadv 

divergence  associated  with  advection 

1/s 

dt 

*thermo/transport  time  step 

3600. s 

dtdyn 

dynamics/transport  time  step  ( A tdyn ) 

dte 

subcycling  time  step  for  elastic  dynamics  (A te) 

s 

dtei 

l/dte,  where  dte  is  the  EVP  subcycling  time  step 

1/s 

dT  mlt 

[see  ice  albedo.F] 

1.  deg 

dumpfile 

*  output  file  for  restart  dump 

dumpfreq 

*  dump  frequency  for  restarts,  y,  m,  or  d 

y 

dumpfreq_n 

*restart  output  frequency 

l 

dragw 

drag  coefficient  for  water  on  ic e*pM, 

0.00536*rhow  kg/m3 

dxt 

width  of  T  cell  (Ax)  through  the  middle 

m 

dxu 

width  of  U  cell  (Ax)  through  the  middle 

m 

dyt 

height  of  T  cell  (Ay)  through  the  middle 

m 

dyu 

height  of  U  cell  (Ay)  through  the  middle 

m 

dvidtd 

ice  volume  tendency  due  to  dynamics/transport 

m/s 

dvidtt 

ice  volume  tendency  due  to  thermodynamics 

m/s 

E 

ecc 

yield  curve  major/minor  axis  ratio,  squared 

4. 

eice(n) 

energy  of  melting  of  ice  per  unit  area  (in  category  n) 

J/rrr 

emissivity 

emissivity  of  snow  and  ice 

0.95 

eps04 

a  small  number 

io  -4 

epsl  1 

a  small  number 

io  -11 

epsl2 

a  small  number 

10  -12 

epsl  3 

a  small  number 

10  -13 

epsl  5 

a  small  number 

10  45 

esno(n) 

energy  of  melting  of  snow  per  unit  area  (in  category 
n) 

J/m2 

evap 

evaporative  water  flux 

kg/nr/s 
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Name 

Description 

Default  Values 

evpdamping 

*if  true,  use  evp  damping  procedure  [17] 

F 

eye 

coefficient  for  calculating  the  parameter  E,  0<eyc<  1 

0.36 

F 

fcor 

Coriolis  parameter 

1/s 

ferrmax 

max  allowed  energy  flux  error  (thermodynamics) 

lxl O'3  W/m2 

fhnet 

net  heat  flux 

W/m2 

fhnet  hist 

net  heat  flux  to  ocean  (fhnet )  for  history 

W/m2 

flat 

latent  heat  flux 

W/m2 

floediam 

effective  flue  diameter  for  lateral  melt 

300.  m 

flw 

incoming  longwave  radiation 

W/m2 

flwout 

outgoing  longwave  radiation 

W/m2 

frain 

rainfall  rate 

kg/m2/s 

frazil 

frazil  ice  growth 

m 

fresh 

fresh  water  flux  to  ocean 

kg/m2/s 

fresh  hist 

fresh  water  flux  (fresh )  for  history 

kg/m2/s 

frzmlt 

freezing/melting  potential 

W/m2 

frz  onset 

day  of  year  that  freezing  begins 

fsalt 

net  salt  flux  to  ocean 

kg/m2/s 

fsalt  hist 

salt  flux  to  ocean  (fsalt )  for  history 

kg/m2/s 

fsens 

sensible  heat  flux 

W/m2 

fsnow 

snowfall  rate 

kg/m2/s 

fsnowrdg 

snow  fraction  that  survives  in  ridging 

0.5 

fsw 

incoming  shortwave  radiation 

W/m2 

fswabs 

absorbed  shortwave  radiation 

W/m2 

fswthru 

shortwave  penetrating  to  ocean 

W/m2 

fswthru  hist 

shortwave  penetrating  to  ocean  (fswthru)  for  history 

W/m2 

fyear 

current  data  year 

fyear  final 

last  data  year 

fyear  init 

*  initial  data  year 

G 
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Name 

Description 

Default  Values 

gravit 

gravitational  acceleration 

9.80616  m/s2 

gridfile 

*input  file  for  grid  info 

grid_type 

* ‘rectangular’  or  ‘displace_pole’  or  ‘column’ 

displaced_pole 

Gstar 

used  to  compute  ridging  participation  function 

0.15 

H 

hfrazilmin 

minimum  thickness  of  new  frazil  ice 

0.05  m 

hi  min 

minimum  ice  thickness  for  thinnest  ice  category 

0.01  m 

hicen 

ice  thickness  in  category  n 

m 

hin  max 

category  limits 

m 

hist_avg 

*if  true,  write  averaged  data  instead  of  snapshots 

T 

histfreq 

*history  output  frequency;  y,  m,  w,  d  or  1 

m 

history  dir 

*path  to  history  output  files 

history  file 

*  output  file  for  history 

hm 

land/boundary  mask,  thickness  (T-cell) 

hmix 

ocean  mixed  layer  depth 

20  m 

hsnomin 

minimum  thickness  for  which  Ts  is  computed 

1.  x  10'6  m 

Hstar 

determines  mean  thickness  of  ridged  ice 

25.  m 

HTE 

length  of  eastern  edge  (Ay)  of  T-cell 

m 

HTN 

length  of  northern  edge  (Ax)  of  T-cell 

m 

HTS 

length  of  southern  edge  (Ax)  of  T-cell 

m 

HTW 

length  of  western  edge  (  Ay)  of  T-cell 

m 

I 

iOvis 

fraction  of  penetrating  visible  solar  radiation 

0.70 

icells 

number  of  grid  cells  with  specified  property  (for 
vectorization) 

ice  ref  salinity 

reference  salinity  for  ice-ocean  exchanges 

4.  psu 

iceruf 

ice  surface  roughness 

5.  x  10’4  m 

icetmask 

ice  extent  mask  (T-cell) 

iceumask 

ice  extent  mask  (U-cell) 

idate 

the  date  at  the  end  of  the  current  time  step 
(yyyymmdd) 

ierr 

general-use  error  flag 
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Name 

Description 

Default  Values 

i(J)hi 

last  i(j)  index  of  physical  domain  (local) 

i(j)l° 

first  i(j)  index  of  physical  domain  (local) 

ilyrl 

index  of  the  top  layer  in  each  cat  (for  eicen ) 

ilym 

index  of  the  bottom  layer  in  each  cat  (for  eicen) 

i(j  )mt_global 

number  of  physical  gridpts  in  x(y)  direction,  local 
domain 

i(j)mt_local 

total  no.  of  gridpoints  in  x(y)  direction,  local  domain 

int  kind 

definition  of  an  integer 

kind(l) 

ipj'P 

local  processor  coordinates  for  writing  debugging 
data 

istep 

local  step  counter  for  time  loop 

istepO 

*number  of  steps  taken  in  previous  run 

0 

istep  1 

total  number  of  steps  at  current  time  step 

K 

kappav 

visible  extinction  coeff.  In  ice,  wavelength  <700 
nm 

1.4/m 

kappan 

visible  extinction  coeff.  In  ice,  wavelength  >700 
nm 

17.6/m 

kcatbound 

*category  boundary  formula 

0 

kdyn 

*type  of  dynamics  (1  =  EVP,  0  =  off) 

1 

kg_to_g 

kg  to  g  conversion  factor 

1000. 

kice 

thermal  conductivity  of  fresh  ice 

2.03  W/m/deg 

kitmin 

minimum  conductivity  of  saline  ice 

W/m/deg 

kitd 

*type  of  ITD  conversions  (1  =  delta  fxn,  1  =  linear 
remap) 

1 

kmt  file 

*  input  file  for  land  mask  info 

krdg_partic 

*ridging  participation  function 

1 

krdgredist 

*ridging  redistribution  function 

1 

ksmooth 

*  1  =  smooth  the  ice  strength 

0 

ksno 

thermal  conductivity  of  snow 

0.30  W/m/deg 

kstrength 

*ice  strength  formulation  (1  =  [35],  0  =  [15]) 

1 

L 

1  conservation  check 

if  true,  check  conservation 
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Name 

Description 

Default  Values 

Lfresh 

latent  heat  of  melting  of  fresh  ice  =  Lsub=Lvap 

J/kg 

lhcoef 

transfer  coefficient  for  latent  heat 

logkind 

definition  of  a  logical  variable 

kind(.true.) 

Lsub 

latent  heat  of  sublimation  for  fresh  water 

2.835  x  106  J/kg 

Lvap 

latent  heat  of  vaporization  for  fresh  water 

2.501  x  106  J/kg 

M 

m  to  cm 

meters  to  cm  conversion 

100. 

ml 

constant  for  lateral  melt  rate 

1.6x1  O'6  m/s  deg"'"2 

m2 

constant  for  lateral  melt  rate 

1.36 

m2  to  km2 

m2  to  km2  conversion 

1  xlO"6 

mask  n(s) 

northern  (southern)  hemisphere  mask 

master  task 

task  ID  for  the  controlling  processor 

mday 

day  of  the  month 

meltb 

basal  ice  melt 

m 

meltl 

lateral  ice  melt 

m 

meltt 

top  ice  melt 

m 

melt  onset 

day  of  year  that  surface  melt  begins 

month 

the  month  number 

MPICOMMICE 

communicator  for  ice  model  internal 
communications  (MPI) 

mpstocmpdy 

m  per  s  to  cm  per  day  conversion 

8.64  x  106 

mpstocompyr 

m  per  s  to  cm  per  year  conversion 

mtask 

local  processor  number  that  writes  debugging  data 

my  task 

task  ID  for  the  current  processor 

N 

nbr  <dir> 

processor  numbers  for  the  N,  S,  E,  W  neighbor 
processors 

neat 

number  of  ice  categories 

5 

ndte 

*number  of  subcycles 

120 

ndyndt 

*number  of  dynamics/advection  steps  under  thermo 

1 

new_day 

flag  for  beginning  new  day 
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Name 

Description 

Default  Values 

new  month 

flag  for  beginning  new  month 

new_week 

flag  for  beginning  new  week 

new  year 

flag  for  beginning  new  year 

ngroups 

number  of  groups  of  flux  triangles  in  remapping 

5 

nilyr 

number  of  ice  layers 

4 

npt 

*total  number  of  time  steps  (dt) 

24 

ntilay 

sum  of  number  of  layers  in  all  categories 

ntracer 

number  of  tracers  transported  in  remapping 

nudiag 

unit  number  for  diagnostics  output  file 

6 

nudump 

unit  number  for  dump  file  for  restarting 

50 

nuforcing 

unit  number  for  forcing  data  file 

49 

nugrid 

unit  number  for  grid  file 

11 

nu  kmt 

unit  number  for  land  mask  file 

12 

nu  nml 

unit  number  for  namelist  input  file 

21 

nu  restart 

unit  number  for  restart  input  file 

50 

nu  rst_pointer 

unit  number  for  pointer  to  latest  restart  file 

52 

numghostcells 

no.  of  rows  of  ghost  cells  surrounding  each 
subdomain 

1 

nyr 

year  number 

o 

oceanmixed  file 

*data  file  containing  ocean  forcing  data 

oceanmixed  ice 

*if  true,  use  internal  ocean  mixed  layer 

T 

ocn  data  dir 

*directory  for  ocean  forcing  data 

omega 

angular  velocity  of  Earth 

7.292  x  1  O'5  rad/s 

one 

array  of  ones  which  is  often  useful 

1. 

opening 

rate  of  ice  opening  due  to  divergence  and  shear 

1/s 

P 

pOOl 

1/1000 

pOl 

1/100 

p027 

1/36 

p055 

1/18 
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Name 

Description 

Default  Values 

Pi 

1/10 

pill 

1/9 

pl5 

15/100 

pl66 

1/6 

P2 

1/5 

p222 

2/9 

p25 

1/4 

p333 

1/3 

p4 

2/5 

p5 

Vi 

p52083 

25/48 

p5625m 

-9/16 

p6 

3/5 

p666 

2/3 

Pi 

TC 

pih 

7r/2 

pi2 

2tc 

pointer  file 

*input  file  for  restarting 

potT 

atmospheric  potential  temperature 

K 

precipunits 

^liquid  precipitation  data  units 

print  global 

*if  true,  print  global  data 

F 

print_points 

*if  true,  print  point  data 

F 

Pstar 

ice  strength  parameter 

2.75xl04  N/m 

puny 

a  small  positive  number 

1  x  10'11 

Q 

Qa 

specific  humidity  at  1 0  m 

kg/kg 

qdp 

deep  ocean  heat  flux 

W/m2 

qqqice 

for  saturated  specific  humidity  over  ice 

1.16378  x  107  kg/m3 

qqqocn 

for  saturated  specific  humidity  over  ocean 

6.275724  x  106 
kg/m3 

Qref 

2  m  atmospheric  reference  specific  humidity 

kg/kg 

R 
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Name 

Description 

Default  Values 

radtodeg 

degree-radian  conversion 

180/n 

radius 

earth  radius 

6.37  x  106  m 

real  kind 

definition  of  single  precision  real 

selected  real  kind(6 
) 

restart 

*if  true,  initialize  using  restart  file  instead  of 
defaults 

T 

restart  dir 

*path  to  restart/dump  files 

restore  sst 

*restore  SST  to  data 

rhoa 

air  density 

kg/m3 

rhofresh 

density  of  fresh  water 

1000.0  kg/m3 

rhoi 

density  of  ice 

917.  kg/m3 

rhos 

density  of  snow 

330.  kg/m3 

rhow 

density  of  seawater 

1 026  kg/m3 

milyr 

real  (nlyr) 

rside 

fraction  of  ice  that  melts  laterally 

s 

saltmax 

max  salinity,  at  ice  base 

3.2  ppm 

sec 

seconds  elapsed  into  idate 

secday 

number  of  seconds  in  a  day 

86400. 

shear 

strain  rate  11  component 

1/s 

shcoef 

transfer  coefficient  for  sensible  heat 

sigl(2) 

principal  stress  components  (diagnostic) 

sinw 

sine  of  the  turning  angle  in  water 

0. 

snoice 

snow-ice  formation 

m 

snowpatch 

length  scale  for  parameterizing  nonuniform  snow 
coverage 

0.02  m 

spval 

special  value  (generally  over  land  or  undefined 
regions,  in  place  of  0) 

1030 

ss_tltx(y) 

sea  surface  slope  in  the  x(y)  direction 

m/m 

sss 

sea  surface  salinity 

psu 

sss_data_type 

*source  of  surface  salinity  data 

sst 

sea  surface  temperature 

C 

sst_data_type 

*source  of  surface  temperature  data 
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Name 

Description 

Default  Values 

stefan-boltzmann 

Stefan-Boltzmann  constant 

5.67  x  10'8  W/nrK4 

stopnow 

if  1 ,  end  program  execution 

strairx(y) 

stress  on  ice  by  air,  in  the  x(y) -direction  (centered  in 

U  cell) 

N/m2 

strairx(y)T 

stress  on  ice  by  air,  x(y)-direction  (centered  in  T 
cell) 

N/  m2 

strength 

ice  strength  (pressure) 

N/m 

stressp 

internal  ice  stress,  an  +  a22 

N/m 

stressm 

internal  ice  stress,  an  -  a22 

N/m 

stress  12 

internal  ice  stress,  ai2 

N/m 

strintx(y) 

divergence  of  internal  ice  stress,  x(y) 

N/m2 

strocnx(y) 

ice-ocean  stress  in  the  x(y)-direction  (U-cell) 

N/m2 

strocnx(y)T 

ice-ocean  stress  in  the  x(y)- dir.  (T-cell) 

N/m2 

strtlx(y) 

surface  stress  due  to  sea  surface  slope 

N/m2 

swv(n)dr(f) 

incoming  shortwave  radiation,  visible  (near  IR), 
direct  (diffuse) 

W/m2 

T 

Tair 

air  temperature  at  1 0  m 

K 

tarea 

area  of  T-cell 

m2 

tarean 

area  of  northern  hemisphere  T-cells 

m2 

tarear 

1 1  tarea 

1/m2 

tareas 

area  of  southern  hemisphere  T-cells 

m2 

If 

freezing  temperature 

C 

Tffresh 

freezing  temp  of  fresh  ice 

273. 15K 

time 

total  elapsed  time 

s 

time  fore 

time  of  last  forcing  update 

s 

Timelt 

melting  temperature  of  ice  top  surface 

0.  C 

tinyarea 

puny  *  tarea 

2 

m 

TLAT 

latitude  of  cell  center 

radians 

TLON 

longitude  of  cell  center 

radians 

tmask 

land/boundary  mask,  thickness  (T-cell) 

tmass 

total  mass  of  ice  and  snow 

kg/m2 
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Name 

Description 

Default  Values 

Tmin 

minimum  allowed  internal  temperature 

i 

O 

O 

O 

Tref 

2m  atmospheric  reference  temperature 

K 

trestore 

*SST  restoring  time  scale 

days 

Tsfc(n) 

temperature  of  ice/snow  top  surface  (in  category  n) 

C 

Tsf  errmax 

max  allowed  Tsjc  error  (thermodynamics) 

5.  x  10'4  deg 

Tsmelt 

melting  temperature  of  snow  top  surface 

0.  C 

TTTice 

for  saturated  specific  humidity  over  ice 

5897.8  K 

TTTocn 

for  saturated  specific  humidity  over  ocean 

5107.4  K 

u 

uarea 

area  of  U-cell 

m2 

uarear 

1  /uarea 

u(v)atm 

wind  velocity,  x(y) 

m/s 

ULAT 

latitude  of  U-cell  centers 

radians 

ULON 

longitude  of  U-cell  centers 

radians 

umask 

land/boundary  mask,  velocity  (U-cell) 

umin 

min  wind  speed  for  turbulent  fluxes 

1.  m/s 

u(v)ocn 

ocean  current,  x(y)  direction 

m/s 

uvel(vvel) 

x(y)-component  of  velocity 

m/s 

uvm 

land/boundary  mask,  velocity  (U-cell) 

y 

vice(n) 

volume  per  unit  area  of  ice  (in  category  n) 

m 

vonkar 

von  Karman  constant 

0.4 

vsno(n) 

volume  per  unit  area  of  snow  (in  category  n) 

m 

W 

week 

week  of  the  year 

wind 

wind  speed 

m/s 

work_gl 

allocatable,  dbl  kind  work  array 

work_g2 

allocatable,  dbl  kind  work  array 

work_gr 

allocatable,  real  kind  work  array 
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Name 

Description 

Default  Values 

write  history 

if  true,  write  history  now 

write_ic 

if  true,  write  initial  conditions  now 

work  11 

(imt  local,  jmt  local)  work  array 

work  12 

(imt  local,  jmt  local)  work  array 

work  a 

(ilo:ihi,  jlo:jhi)  work  array 

workb 

(ilo:ihi,  jlo:jhi)  work  array 

write  restart 

if  1,  write  restart  now 

Y 

ycycle 

*number  of  years  in  forcing  data  cycle 

yday 

day  of  the  year 

year  init 

*the  initial  year 

Z 

zlvl 

atmospheric  level  height 

m 

zref 

reference  height  for  stability 

10.  m 

zTrf 

reference  height  for  Tre/,  Qref 

2.  m 

zvir 

gas  constant  (water  vapor)/gas  constant  (air)  -1 

0.606 
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7.0  NOTES 

7.1  Acronyms  and  Abbreviations 


Acronym 

Definition 

CCSM 

Community  Climate  System  Model 

CICE 

Los  Alamos  Sea-Ice  Model 

DMSP 

Defense  Meteorological  Satellite  Program 

DTG 

Date-Time-Group 

EVP 

Elastic-Viscous-Plastic 

FNMOC 

Fleet  Numerical  Meteorology  and  Oceanography  Center 

HYCOM 

HYbrid  Coordinate  Ocean  Model 

I/O 

Input/Output 

ITD 

Ice  Thickness  Distribution  model 

LANL 

Los  Alamos  National  Laboratory 

MPDATA 

Multidimensional  Positive  Definite  Advection  Transport  Algorithm 

MPI 

Message  Passing  Interface 

NAVOCEANO 

Naval  Oceanographic  Office 

NCAR 

National  Center  for  Atmospheric  Research 

NCOM 

Navy  Coastal  Ocean  Model 

NOGAPS 

Navy  Operational  Global  Atmospheric  Prediction  System 

NRL 

Navy  Research  Laboratory 

PIPS 

Polar  Ice  Prediction  System 

POP 

Parallel  Ocean  Program  model 

PSI 

Planning  Systems,  Incorporated 

S 

Salinity 
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SDD 

Software  Design  Description 

SGI 

Silicon  Graphics  Incorporated 

SHEBA 

Surface  Heat  Budget  of  the  Arctic  Ocean 

SSC 

Stennis  Space  Center 

SSM/I 

Special  Sensor  Microwave/Imager 

SSS 

Sea  Surface  Salinity 

SST 

Sea  Surface  Temperature 

SUM 

Software  Users  Manual 

T 

Temperature 
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Appendix  A 


Table  of  Namelist  Options 


Name 

Type/Options 

Description 

Default  Values  / 
Directory  Location 

albicei 

0  <  a  <  1 

near  infrared  ice  albedo  for 
thicker  ice 

albicev 

0  <  a  <  1 

visible  ice  albedo  for  thicker 
ice 

albsnowi 

0  <  a  <  1 

near  infrared,  cold  snow 
albedo 

albsnowv 

0  <  a  <  1 

visible,  cold  snow  albedo 

advection 

remap 

linear  remapping  advection 

’remap’ 

mpdata 

2nd  order  MPDATA 

upwind 

1st  order  MPDATA 

atm  data  dir 

path/ 

path  to  atmospheric  forcing 
data  directory 

atm  data  type 

default 

constant  values  defined  in 
the  code 

near 

NCAR  bulk  forcing  data 

LYq 

AOMIP/Large- Y  eager 
forcing  data 

dbug 

true/false 

if  true,  write  atm/ocn  data 
diagnostics 

.false. 

diagfile 

filename 

diagnostic  output  file 

diag_type 

stdout 

write  diagnostic  output  to 
stdout 

‘stdout’  (if  uncoupled) 

file 

write  diagnostic  output  to 
file 

diagfreq 

integer 

frequency  of  diagnostic 
output  in  dt 

30 

eg-,  10 

once  every  10  time  steps 

dt 

seconds 

thenno/transport  time  step 
length 

2800. 

dumpfile 

filename  prefix 

output  file  for  restart  dump 

’iced’ 

dumpfreq 

y,  m,  d 

write  restart  every 
dumpfreq _n  for  years, 

'd' 
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Name 

Type/Options 

Description 

Default  Values  / 
Directory  Location 

months,  days 

dumpfreq_n 

integer 

frequency  restart  data  is 
written 

evpdamping 

true/false 

if  true,  damp  elastic  waves 
[61 

.false. 

fyear  init 

yyyy 

first  year  of  atmospheric 
forcing  data 

f  <var> 

true/false 

write  <var>  to  history 

grid_fde 

filename 

name  of  grid  file  to  be  read 

‘grid’ 

grid_type 

rectangular, 

displaced_pole 

rectangular:  defined  in 
rectgrid 

displaced_pole:  read  from 
file  in  popgrid 
pipsgrid:  read  from 
file  in  pipsgrid 

‘pipsgrid’ 

hist_avg 

true/false 

write  time-averaged  data  if 
true 

write  snapshots  of  data  if 
false 

.false. 

hist  dir 

path/ 

path  to  history  output 
directory 

histfreq 

y,  m,  w,  d,  1 

write  history  output  once  a 
year,  month,  week,  day,  or 
every  time  step 

V 

history  file 

filename  prefix 

output  file  for  history 

'iceh' 

ice_ic 

default 

latitude  and  SST  dependent 

‘default’ 

none 

no  ice 

istepO 

integer 

initial  time  step  number 

0 

kcatbound 

0/1 

if  0,  original  category 
boundary  formula 
if  1 ,  new  category  boundary 
formula 

1 

kdyn 

0/1 

if  0,  EVP  dynamics  OFF 
if  1 ,  EVP  dynamics  ON 

1 

kitd 

0/1 

if  0,  delta  function  ITD 
approx. 

if  1 ,  linear  remapping  ITD 
approx. 

1 

kmt  file 

filename 

name  of  land  mask  file  to  be 

'kmt' 
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Name 

Type/Options 

Description 

Default  Values  / 
Directory  Location 

read 

krdg_partic 

0/1 

if  0,  old  ridging  participation 
function 

if  1 ,  new  ridging 
participation  function 

1 

krdgredist 

0/1 

if  0,  old  ridging 
redistribution  function 
if  1 ,  new  ridging 
redistribution  function 

0 

kstrength 

0/1 

if  0,  [15]  ice  strength 

formulation 

if  1,  [35]  ice  strength 

formulation 

1 

ndte 

integer 

number  of  EVP  subcycles 

120 

ndyndt 

integer 

number  of 

dynamics/ advection/  ridging 
steps  per  thermo  timestep 

1 

npt 

integer 

total  number  of  time  steps  to 
take 

oceanmixed  file 

filename 

data  file  containing  ocean 
forcing  data 

oceanmixed  ice 

true/false 

active  ocean  mixed  layer 
calculation 

.true,  (if  uncoupled) 

ocn  data  dir 

path/ 

path  to  oceanic  forcing  data 
directory 

'/scr/posey/pips3c/data_in/' 

print_points 

true/false 

print  diagnostic  data  for  two 
grid  points 

.true. 

precipunits 

mm_per  month 

liquid  precipitation  data 
units 

mm_per  sec 

(default;  MKS  units) 

print  global 

true/false 

print  diagnostic  data,  global 
sums 

.true. 

pointer  file 

pointer  filename 

contains  restart  filename 

restart 

true/false 

initialize  using  restart  file 

.true. 

restart  dir 

path/ 

path  to  restart  directory 

restore  sst 

true/false 

restore  SST  to  data 

sss_data_type 

default 

constant  values  defined  in 
the  code 

clim 

climatological  data 
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Name 

Type/Options 

Description 

Default  Values  / 
Directory  Location 

near 

POP  ocean  forcing  data 

sst_data_type 

default 

constant  values  defined  in 
the  code 

clim 

climatological  data 

near 

POP  ocean  forcing  data 

trestore 

integer 

SST  restoring  time  scale 
(days) 

ycycle 

integer 

no.  of  years  in  forcing  data 
cycle 

year  init 

yyyy 

the  initial  year,  if  not  using 
restart 
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